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THE PROGENY OF AN ENTIRE POPULATION! 


By ALFrep J. LOTKA 


Metropolitan Life Insurance Company 


The literature on renewal theory has grown to considerable dimensions, 
until even admittedly incomplete bibliographies list over 100 titles. But a 
surprisingly small proportion of these publications exhibits any practical ap- 
plications to concrete data, and such applications as have been made (e.g. by 
Wicksell, Hadwiger, Rhodes) are for the most part of restricted scope. 

Anyone who has been following the development will, I think, feel that this is 
unfortunate. It has a double disadvantage. On the one hand the purely 
theoretical discussions emphasize difficulties which in practice may be relatively 
unimportant, being inherent either in some of the unrealistic ad hoc examples 
discussed, or in the expressions used to fit smooth curves to the basic data, 
rather than in these data themselves. On the other hand some real difficulties 
in application to actual data seem to require further clarification. 

Several of the applications that have been made, including some of my own, 
are restricted to following up the “progeny”’ of a “‘population element’’ com- 
prising only individuals all originating at the same time and therefore all of the 
same age (in the case of industrial equipment installation all made at one point 
of time). The analysis set forth in the treatment of this special case is competent 
also to deal with the practically more important case of the progeny of an initial 
population of given age distribution, though no example of this has hitherto 
been published.? Such an example will now be given, and at the same time this 
will afford an opportunity to clarify some points in the presentation of the more 
general case. 

Let N; be the total number of females at time ¢, and c,(a) the number com- 
prised within the age limits a and a + da. Also, let m,(a) be the age-specific 
fertility of females of age a, counting daughters only. If a and w are, respec- 
tively the lower and the upper limit of the female reproductive period, and B(t) 
the annual births of females, then 


(1) Bit) = | : N,¢,(a)m,(a) da. 


However, it is not in this perfectly general form that the relation is to be ap- 
plied. The case to be considered is that in which the “initial” population is 
throughout its ‘future’ development, subject to constant age-specific fertility 


1 Compare A. J. Lotka, ‘‘The progeny of a population element,’’ Am. Jour. Hygiene, 
Vol. 8 (1928), p. 875. 

2 An example was given by the writer in an oral communication to the Eighth American 
Scientific Congress, May 1940, the Proceedings of which have not so far been published. 
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and mortality. If we denote the “initial” time by ¢ = w (which we can do since 
the zero of time is arbitrary), we can then write 


(2) Bit) = f N,c,(a)m,(a) da, t> wo. 


Also, if p:-.(a) is the probability for a female born at time 7 = ¢t — a of sur- 
viving to time ¢t, being then a years old, we have 


(3) B(t — a)p:.(a) = N,c,(a), 


and, in particular, since in the case under consideration p,;_,(a) is constant for 
t — a > ao, i.e., for individuals born after t = w 


(4) B(t — a)p.(a) = N,c.(a), t>a+t+u. 


Now, we have been at liberty for the ‘‘future” values of m,(a) and p;_,(a) 
to make the arbitrary assumption that they retain their values as of t = w and 
t — a > w, respectively. But for the “past” of the system under consideration 
we do not have equal liberty, for any assumption we make must be compatible 
with 

(a) the initial age distribution 
(b) equation (1). 

We can, however, within these limitations, assume that (4) still holds for 
0 <t < w, thus 
(5) B(t — a)p.(a) = Nrex(a), t>0. 


Introducing this in (1) we have 


(6) Bit) = [ B(t — a)p.(a)m,(a) da, t > 0. 


But we cannot now, further assume that 
(7) m,(a) = m,(a), t > 0, 
for, in general, this would make (6) incompatible with (1). 

We can, however, split the integral in (6) into two parts, thus 


(8) Bi) = | " Be ~ daldnidds + [ Bt — dadanded de, 


with the assumption, only in the range a <t, 

(9) . m,(a) = m,(a), a<t. 
Denoting the first integral in (8) by F(é), and contracting p.(a)m,(a) to 

¢.(a), we may write (8) in the form 


(10) B(t) 


FW) + [Be = ae.(a) da, 
F() + 6), 


(11) 
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with 

F(t) =0 t>w 
(12) 

\F(t) = Bit) 0<t<a 
and 
(13) | Bi) = | 8 - dalde& t>w! 


The assumption (9) has a definite physical meaning. The integral in (6) 
has been so split that the first part, F(t), gives the births of daughters from 
mothers who themselves were born before t = 0, while the second part, A(t), 
gives the births of daughters from mothers born after ¢ = 0. Equation (9) 
therefore expresses the assumption that for mothers born at or after t = 0, 
the age-specific fertilities for ages a < t have the same values m,(a), independent 
of ¢, as prevail for = w. But at time ¢ there are no mothers of age a > t, 
who were born after ¢ = 0. Hence the assumption (9) can be quite simply 
stated to the effect that the age-specific fertilities m.(a) apply to all mothers 
born after time = 0. This assumption cannot, in general be made for mothers 
born before t = 0, because it would not, in general, be compatible with the 
given initial age distribution and at the same time with assumption (5). Hence 
in the first integral of (8), denoted by F(t) in (10), we must write m,(a), not 
m,(a). 

Equation (10) is of the form discussed by G. Herglotz,* who writes its solution, 
for t > 0, in the form of an exponential series. 


(14) B(t) = 2Q;e** 
where the exponents 7; are the roots of the characteristic equation, 
(15) ar) = [ yl) da = 1, 


while the coefficients Q; are given by 


I ” Fe! dt 
(16) Q, = ———___—__. 
/ ae’ (a) da 


There is only one real root of (14), since g(a) = 0, for all values of a. For 
complex roots it is convenient to write the corresponding terms of the series (14) 
in trigonometric form 


(17) Qe = 2Ue"' cos vt — 2Ve" sin vt, 
(18) 2./(U? + V2%)e"* cos (vt + 8), 


3 Since g(a) = Ofora> w. 
4 Math. Annalen, Vol. 65 (1908), pp. 87 et seq. 
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where 
tan @ = V/U, 


U 
(19) cos § = Jt? + v2? 


s 


sin @ = VU dv 


and 










- RG+SH 
- ce 
- RH—-SG 
2) - 
in which 
(22) G = [ ae “° cos va ¢,(a) da, 
(23) H = [ ae “* sin va ¢,,(a) da 
(24) r=[ e “' cos vt F(t) dt, 
0 
(25) S= : e “' sin vt F(t) dl. 
0 


For purposes of numerical application to the problem here censidered, we must 
express the annual births B(t) for ¢ < w in terms of the given “initial” age 
distribution at time w. 


We have, generally 


/ 


N,c,(a) Nucla + w — 
26 BG — a) = —— = ——— =» 
- oO Me” Metro 8 
since individuals of age a at time t, area + w — ¢t years old at time w. 
Introducing the relation (26) in (10) we have 












(27) BO = FO+ f “pe Fao Plaiele de 





and 
a - Nwcu(a + w — 
(28) F() = BY [ Nett FO pula) da, 
a Culw — t) | ca(a + w — t) 
(29) — ee [* Pala > aon t) (a) da, 


(11a) = Bit) — BO). 
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Note that, in computing the integral 8(t) for any particular value of ¢, the 
argument of the function c, runs from a +w— ¢ttow. Thus, for example, if 
the zero of time is 1865 and ¢ = w is at 1920, then, in computing F(35), i.e., 
the value of F for 1900, the range of the argument of c, in the integral will be 
from 10 + 55 — 35 to 55, i.e., from 30 to 55. 

Numerical Example. By way of a numerical illustration these principles will 
now be applied to a concrete case. We shall start with the age distribution of 
the white female population of the United States as constituted in 1920, for 
which previous publications furnish some of the required data, including the real 
root and the first three pairs of complex roots of the characteristic equation. 

From this “‘initial’’ age distribution in 1920 it is necessary first of all to com- 
pute the auxiliary function F(t) for the 55 vears prior to 1920. The first term 
B(t) in the right hand member of (28) is very easily computed for successive 
values of ¢ from the relation (5a), which simply expresses the fact that persons 
a vears old in the year w, i.e., 1920, are the survivors of the B(w — a) persons 
born in the year w — a. 

(5a) N.c.(a) = Bw — a)p.(a). 

In the diagram Fig. 1, which is drawn in stereographic projection, the age 
distribution of the (white female) population of the United States in 1920 is 
represented as plotted in a plane reaching forward at right angles to the plane 
of the paper. Successive points of B(t) for 0 < t < w, have been computed 
“by survivals’ according to (5a) and plotted as a curve in the plane of the 
paper “‘at the back” of the diagram. The arrows indicate for a selected point, 
namely age 25 in 1920, the path of the computation according to equation (5a.) 

The second term @(¢) in the expression (lla) for F(t) was computed from the 
age distribution in 1920, the rates of survival from previous years into 1920,° 
and the age-specific fertility at each age in the reproductive period, 10 to 55, 
on the basis of the relation (28). The results, for this second term in the ex- 
pression for F(t) computed for every fifth calendar year back of 1920 to 1875 
and interpolated for intervening years,’ were also plotted as a curve in the rear 
plane of the diagram. The shaded area in the curve for the age distribution in 
1920, and the arrows leading from this shaded area to the curve 


(10, 11) a) = | BO - ae.(a) da 
(29, 11a) — [ Neca + w — 2) (a) da, 


a Pu(a Sa = t) ~ 
indicate in this case the path of the computation according to, equation (28). 


5 Using the Foudray life table for white females in 1919-1920. In the first quinquennial 
age group, the following values were used: 
p(0.5) = .9460 p(2.5) = .9135 
p(1.5) = 9235. p(3.5) = .9080 p(4.5) = .9040 


6 This term vanishes for t < 10, i.c., back of 1875. 
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From these two curves, taking differences, the curve of F(t) = B(t) — B(d) 
was plotted, as shown. 

With the values of F(#) thus obtained, we may proceed, by formulae (14) to 
(25), to compute values of B(t) for all values of t > 0. So far as the period 
1865 to 1920, corresponding to 0 < t < w, is concerned, this merely means that 
we have an analytical expression to fit what is essentially a fundamental datum 
of the problem. For values of ¢ > w the formula gives us a continuation of the 
function B(t) for all future time so long as the given age-specific fertility and 
mortality holds. 


CALENDAR YEAR 
168s 1895 


1905 19 
es B(t)= births computed from survivals ~~] 


MILLIONS 


=births computed from a7 
ity, exclusive of those \ 
whose mothers were born ~\ 
Lippy? Ty 


& 
x 
8 
X 
x 


Fic. 1. Graph illustrating computation of auxiliary function F(t) from “‘initial’’ age 
distribution. 


The final results of this computation are exhibited in Figs. 2, 3 and 4. Of 
these, Fig. 2 exhibits the first, second and third oscillatory components for the 
period from 1890 forward. It will be seen that the waves are heavily damped, 
so that after a relatively short period the aperiodic component dominates the 
course of events. 

Fig. 3 exhibits, for the years from 1865 to 1920, i.e., for the period 0 < t < w, 
the aperiodic component (in a dashed line) and, as indicated by small circles, 
the sum of this component plus the three oscillatory components. It will be 
seen that from about 1890 forward the points so obtained follow rather closely 
the value B(t) derived by survivals from the age distribution in 1920. 
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DEVIATIONS FROM APERICDIC COMPONENTS 


Fic. 2. First three oscillatory components of total annual births 


CALENDAR YEAR 
1865 1875 1685 1895 1905 1915 1920 


MILLIONS 


4 /B(t)=births computed from . 
fertility, exclusive of those 
whose mothers Were born A 


@ Sum of first four components 
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Fic. 3. Graph of functions B(t), B(t), and F(t) for 0 < t < w, i.e., for 1865 to 1920, to- 
gether with aperiodic component; also, summation of aperiodic and first three os- 


cillatory components. 
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Oscillatory Components 
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/ ----Long Time Trend 
(Aperiodic Component) 


© Sum of Aperiodic 
and First Three 
Oscillatory Components 


From |890 up to 1920: 
Births Computed from 
Survivors in 1920 

From (920 forward: 
Births “Projected” on basis of 

-Schedules of Fertility 

and Mortality in 1920 


CALENDAR YEAR 


TABLE I 


Constants of the Series Solution (14) of Integral Equation (10) to Third Oscillatory 
Component Inclusive t = 0 at 1865 
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Second 
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Third 


48.849 X 107 
37 .008 
17.266 
11.684 X 10° 


i— 16.543 X 108 


88.014 X 10° 
48.808 X 10! 






Fic. 4. Sum of aperiodic and three oscillatory terms of series solution compared with 
results of ‘‘step by step’’ computation of annual births. 





—9.804 X 10- 


2 
2 
2 


PROGENY OF A POPULATION 123 


Prior to about 1890, four components alone are quite inadequate, and the 
corresponding points have been omitted from the diagram. The lack of con- 
cordance, with such limited components, is inconsequential in this part of the 
series, since the purpose of this part of the work was merely to compute the 
auxiliary function F(t), and the fit obtained for B(t) in this range, so far as it goes, 
is merely a by-product, the main interest being in the course of B(t) for t > w, 
i.e., in the years following 1920. 

This course is charted in Fig. 4, in which the points obtained by the series 
solution (14) of (10) are again shown as small circles, while the fully drawn curve 
is derived from my previous publication ‘The Progressive Adjustment of Age 
Distribution to Fecundity.”” The annual births in that case were obtained 
“step by step’ by computing age distributions by survivals for successive 


TABLE II 


United States White Female Population 1920, Observed; Also, the Same Projected 
Forward for Later Years® 


Year Population, Births, | Birth rate per 1,000 
thousands thousands per annum 


1920 49 ,390 1,082 23 .32 
1930 51 ,727 1,162 22.46 
1940 56 ,910 1,252 22.00 
1950 61 ,639 1,307 21.20 
1960 65 ,835 1,379 20.95 
1970 69 ,829 1 ,465 20.98 
1975 71,828 1,504 20.94 
1980 | 73 ,850 1,543 20.89 
1985 | 75 ,902 1,584 | 20.87 











quinquennial periods, and applying to the reproductive age groups, in each 
case, the values of the reproductivity m,(a). 

It will be seen that the points obtained by the solution (14) follow very closely 
those computed “step by step,” although in the computation of the latter an 
approximation was made, using pivotal values of p.(a) for the several quin- 
quennial age groups. A slight error introduced in this way would tend to be 
cumulative, and perhaps accounts for the fact that towards the end of the 
period covered (1985), the two sets of values diverge slightly. Even so, in 1985, 
the divergence is only about .4 percent. 

The series solution has, of course, the advantage that it gives directly the 
result for any particular point of time, whereas the “step by step” method re- 


7 Jour. Washington Acad. Sci., Vol. 16 (1926), p. 505. 
8 Calculated step by step from survival ratios and age specific fertilities, both held 


constant as of 1920 (reproduced for ready reference from Jour. Wash. Acad. Sci., Vol. 16, 
p. 505). 
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quires the computation of the annual births for all intervening points in order 
to obtain the result for the chosen point of time. 

Furthermore, the series tells us at once that the course of events is of the nature 
of a trend proceeding in geometric progression upon which are superposed a 
series of damped oscillations, of which the fundamental has a wave length equal 
approximately to the mean length of one generation from mother to daughter, 
i.e., about 28 years. 

Alternative procedure. The procedure set forth in the preceding sections in- 
volves not only arbitrary assumptions regarding the values of p(a) and m(a) 
for “future” time, which are fundamental to the problem under consideration, 
but involves further incidental assumptions regarding their values prior to the 
“initial” condition at the instant denoted by t = w. These incidental assump- 
tions are in a sense superfluous, since the future history of the system is com- 
pletely determined by the initial age distribution and the assumed “future” 
values of p(a) and m(a). The additional assumptions were introduced merely 
for the purpose of translating the initial age distribution into a series of values of 
B(t) for 0 < t < w,i.e., prior to the given initial age distribution. 

In actual fact the age distribution at time ¢ = w did not arise in the manner 
assumed; actually both p(a) and m(a) undoubtedly varied in the period 1865 
to 1920, and migration also affected the situation. The quantity F(t) intro- 
duced in equation (10) is, in fact, a purely auxiliary function having no direct 
relation to the biological events at time t < w. 

An alternative procedure which would avoid these conflicts, and introduce 
assumptions only regarding “future”? values of p(a) and m/(a), would be to 
compute B(t) step by step over the period from B(1920) to B(1920 + w) = 
B(1975). 

Placing the zero of time ¢ = 1920 this would give B(t) forO0 <t<w. For 
t > w we should have, simply 


Bit) = / B(t — a)¢is20(a) da, t >, 
using, in the evaluation of the integral, the values of B(t — a) obtained by the 


step by step process. 
We could here also split the integral into two parts 


BW) = [Be — adewnla) da + [Bit — a)ewaa) da 


=F) + [BU aeun(a) da. 


But the function ¢920(a) is now the same in the two integrals, and there is no 
occasion, in this case, for distinguishing the two parts of the integral. 
If this procedure is adopted, its application to the course of B(t) for t > w, 
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i.e. beyond 1975, is of minor interest, for by that time it has practically settled 
down to the aperiodic (exponential) component, the oscillations being greatly 
damped down. The major interest in the result of a computation carried out 
by this procedure would be in the fitting of a series of the form (14) to the 
function B(t) in the range 1920 to 1975, which, in this setting, figures as a known 
“arbitrary” function. 


Of the two alternative procedures the one carried out in detail in the text 
and the numerical example is of greater interest, as exhibiting in greater gene- 
rality the application of the Hertz-Herglotz solution. 


BIBLIOGRAPHY 


A few more titles may be added to previously published bibliographies (see these Annals, 
Vol. 10 (1939), p. 22; Vol. 12 (1941), p. 266). 


1935 W. Méscuter, “Untersuchungen iiber Eintrittsgewinn und Fehlbetrag einer Ver- 
sicherungskasse,’’ Bull. de l’Association des Actuaires Suisses, October (1935), 
p. 129. 

1936 ALi AFzaLipour, Contribution a l’étude de la théorie mathématique de la démographie, 
Thése, 1936. 

1938 H. Hapwicer, ‘‘Ein Konvergenzkriterium fiir Erneuerungszahlen,’’ Skandinavisk 
Aktuarietidskrift, Vol. 3-4, p. 226. 

H. Hess, Anwendung der logistischen Funktion in der mathematischen Bevilkerungs- 
theorie, Inaugural Dissertation der philosophischen Fakultat der Universitat 
Bern. 

W. M. Dawson and W. Epwarps Demine, ‘“‘On the problem of natural increase,” 
Growth, Vol. 2, p. 319. 

A. W. Brown, ‘‘A note on the use of a Pearson Type III function in renewal theory,” 
Annals of Math. Stat., Vol. 11 (1940), p. 448. 

E. Zwineoat, ‘‘Entwicklung von Personengesamtheiten—Zusammenfassender Be- 
richt,’’ Twelfth International Congress of Actuaries, Vol. 3, p. 263. 

E. KEINnaNnEN, ‘Uber die Altersverteilung der Bevélkerung,” Twelfth International 
Congress of Actuaries, Vol. 3, p. 305. 

1940 S. S. Townsenp, ‘‘Some observations on the internal variation in groups of lives 
assured in industrial assurance and in the general population of England and 
Wales,”’ Twelfth International Congress of Actuaries, Vol. 3, p. 319. 

1940 R. Tarsan, ‘Untersuchungen tiber den Kapitalbedarf des Lebensversicherungs- 
geschiftes,’’ Twelfth International Congress of Actuaries, Vol. 3, p. 335. 

1940 M. Pressurcer, ‘Sur l’étude générale des collectivités de personnes,” Twelfth 
International Congress of Actuaries, Vol. 3, p. 353. 

1940 A. Maret, ‘“‘Direkte Berechnung der Vorgangsfunktionen einer offenen Gesamtheit,’’ 
Twelfth International Congress of Actuaries, Vol. 3, p. 387. 

1940 E. Zwinecr, “Uber Zusammenhinge zwischen der technischen Stabilitat einer 
Sozialversicherungskasse und der Entwicklungsformel fiir den Versicherten- 
bestand,” Twelfth International Congress of Actuaries, Vol. 3, p. 395. 

1940 H. Hapwicer and W. WecmMUt ter, ‘‘Entwicklung und Umsichichtung von Personen- 
gesamtheiten,”’ Twelfth International Congress of Actuaries, p. 369. 

1940 W. Dossernack and G. T1e7Tz, ‘‘Die Entwicklung von Personengesamtheiten vom 
Standpunkt der Sozialersicherungstechnik,”’ ‘‘Zwélfter Internationaler Kongress 
Der Versicherungs-Mathematiker,’’ Vol. 4, p. 233. 

1941 R.C. Geary, “Irish population prospects considered from the viewpoint of reproduc- 
tion rates,’’ Statistical and Social Inquiry Society of Ireland, 1941. 





126 ALFRED J. LOTKA 


1941 H.Hapwicer, ‘‘Eine Formel der mathematischen Bevélkerungstheorie,’’ Mitteilungen 
der Vereinigung Schweizerischer Versicherungsmathematiker, Vol. 41, p. 67. 

IL. Freraup, ‘‘Le renouvellement, quelques problémes connexes et les equations 
intégrales de cycle fermé,’’ Mitteilungen der Vereinigung Schweizerischer Ver- 
sicherungsmathematiker, Vol. 41, p. 81. 

W. FE.ier, ‘“‘On the integral equation of the renewal theory,’’ Annals of Math. 
Stat., Vol. 12 (1941), p. 243. 

Harro BERNARDELLI, ‘Population waves,” Journal of the Burma Research Society, 
Vol. 31. 





ASYMPTOTICALLY SHORTEST CONFIDENCE INTERVALS! 


By ABRAHAM WALD? 
Columbia University 


The theory of confidence intervals, based on the classical theory of proba- 
bility, has been treated by J. Neyman.’ While Neyman considers the case of 
small samples, we shall deal here with the limit properties of the confidence 
intervals if the number of observations approaches infinity. 


1. Definitions. We will start with some of Neyman’s definitions. Let 
f(x, 6) be the probability density function of a variate z involving an unknown 
parameter @. Denote by E, a point of the n-dimensional sample space of n 
independent observations on x. If p(E,) denotes for each E, a subset of the 
real axis, the symbol P[p(Z,,)c@’ | 6’’] will denote the probability that p(E,) con- 
tains 6’ under the hypothesis that 6’’ is the true value of the parameter. Let 
6(E,,.) and 6(E,) be two real functions defined over the whole sample space such 
that 0(E,) < 6(E,). The interval 5(E,) = [6(E, , 6(E,)] is called a confidence 
interval of 6 corresponding to the confidence coefficient a (0 < a < 1) if 
P[s(E,,)c6 | 6] = a for all values of @. 

The interval function 6(E,,) is called a shortest confidence interval of 6 corre- 
sponding to the confidence coefficient a if 
(a) P[5(E,)cé | 6] = a for all values of 6, and 
(b) for any interval function 6’(Z,) which satisfies the condition (a) we have 

P[5(E,)c@’ | 6”) < P[s'(E,)c@’ | 6”), 
for arbitrary values 6’ and 6’. 

The interval function 6(Z,) is called a shortest unbiased confidence interval 
of 6 if the following three conditions are fulfilled: 
(a) P[é(E,)cé | 6] = ao for all values of @. 

(b) P[5(E,)cé’ | 6’) < a@ for all values of 6’ and @’’. 


(c) For any interval function 5’(EZ,) for which the conditions (a) and (b) are 
satisfied, we have 


P[5(E,)c’ | 6] < Pls’(E,)cé’ | 6”), 
for all values of 6’ and 6”. 


For any relation R we shall denote by P(R | 6) the probability that R holds 
under the hypothesis that @ is the true value of the parameter. Similarly for 


1 Presented at a joint meeting of the Institute of Mathematical Statistics and the Ameri- 
can Mathematical Society in Hanover, September, 1940. 

2 Research under a grant-in-aid from the Carnegie Corporation of New York. 

3 J. NeyMAN, ‘“‘Outline of a theory of statistical estimation based on the classical theory 
of probability,” Phil. Trans. Roy. Soc. London, Vol. 236 (1937), pp. 333-380. 
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any region Q,, of the n-dimensional sample space the symbol P(Q,, | @) will denote 
the probability that the sample point £, falls in Q, under the hypothesis that 6 
is the true value of the parameter. 

In all that follows we shall denote a region of the n-dimensional sample space 
by a capital letter with the subscript n. 

A real function 6(E,) is called a best upper estimate of 6 if the following two 
conditions are fulfilled: 
(a) Pe < 6(E,) | 6] = a@ for all values of 0. 
(b) For any function 6’(E,) which satisfies the condition (a) we have 


Plo’ < 6(E,) | 6] < Pie’ < 6(Z,) | 6”) 


for all values 6’ and @” for which @’ > @”’. 

A real function 9(EZ,) is called a best lower estimate of 6 if the following two 

conditions are fulfilled: 

(a) P(@ > O(E,) | 6] = a for all values of 0. 

(b) For any function 6’(EZ,) which satisfies the condition (a) we have 
P(@’ > OE.) | 0] < Pla’ > o'(E,) | 6") 

for all values of 6’ and @” for which @’ < @’’. 

We will extend the above definitions of Neyman to the limit case when n 
approaches infinity. 

DeFIniTIon I: A sequence of interval functions {5,(E,)} (n = 1, 2,---) ts 
called an asymptotically shortest confidence interval of 6 if the following two conditions 
are fulfilled: 

(a) P[5,(E,)cé | 6] = a for all values of @. 
(b) For any sequence of interval functions {5,(E,)} (n = 1, 2,--+,ad inf.) 
which satisfies (a), the least upper bound of 
P[5,(E,)c0’ | 6] — P[5,(Eq)c6’ | 6”) 
with respect to 6’ and 6” converges to zeroasn — ~. 

Derinition II: A sequence of interval functions {5,(E,)} is called an asymp- 
totically shortest unbiased confidence interval of @ if the following three conditions 
are fulfilled: 

(a) P[é,(E,)cé | 6] = a@ for all values of 0. 
(b) The least upper bound of P{é,(E,)c6’ | 6’’] with respect to 0’ and 6” converges 
toa withn— «. 
(c) For any sequence of interval functions {6,(E,)} which satisfies the conditions 
(a) and (b), the least upper bound of 
P[5,(E,)c’ | 6] — P[d,(En)c@’ | 6”) 
with respect to 6’ and 6’' converges to zero withn > ~. 

Derinition III: A sequence of real functions {6,(E,)} (n = 1, 2, --- , ad inf.) 
is called an asymptotically best upper estimate of 0 if the following two conditions 
are fulfilled: 

(a) P{é < 6,(E,) | 6] = a@ for all values of 0. 
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(b) For any sequence of functions {6,(E,)} which satisfies (a) the least upper 
bound of 
Pe’ < 6,(E,) | 6] — Pl@’ < 6,(E,) | 6”) 
in the domain 6’ > 6” converges to zero withn > ~. 
Derinition IV: A sequence of real functions {6,(E,)} is called an asympto- 
tically best lower estimate of 6 if the following two conditions are fulfilled: 
(a) P(@ > 8,(E.) | 6] = @ for all values of 06. 
(b) For any sequence of functions {6,(E,)} which satisfies (a) the least upper 
bound of 
P[6’ > 8,(En) | 0] — Pl6’ > 9,(Es) | 0” 
in the domain & < 6’’*converges to zero with n — o. 


2. Two Propositions. Proposition I: Let {W,(6)} (n = 1, 2, ---, ad inf.) 
be for each 6 a sequence of regions such that the following two conditions are fulfilled: 
(a) P[W,(6) | 6] = 1 — a for all values of 6. 

(b) For any sequence of regions {Z,(0)} which satisfies (a) the least upper bound of 
P[Z,(6’) | 6”) — P[W,(6’) | 6”) 
in the domain 6’ > 6'(6’ < 6”) converges to zero withn > ~. 


Denote by pn»(E,) the set of all values of 6 for which E,, does not lie in W,(0). Then 
we have 


(c) Plpna(E.)cé | 6] = a for all values of 0. 
(d) For any sequence of set functions {p,(E,)} which satisfies (c), the least upper 

bound of 

Plon(En)c6’ | 6] — Plon(En)cé’ | 6’) 
in the domain 6’ > 6''(6’ < 6) converges to zero withn — ~. 
Proposition II: Let {W,(6)} be for each 0 a sequence of regions such that the 

following three conditions are fulfilled: 
(a) P(W,(@) | 6] = 1 — a for all values of 6. 

(b) The greatest lower bound of P[W,,(6’) | 6’’] converges to1 — awithn > ~. 
(c) For any sequence {W7,(6)} which satisfies (a) and (b), the least upper bound of 
P(W,.(6’) | 6] — P[W,(6’) | 0”) 

with respect to 6’ and 6’ converges to 0 withn > ~. 


Denote by p,(Ex) the set of all values of @ for which E,, does not liein W,(8). Then 
we have 


(d) Plpn(E,)cé | 6] = a for all values of 8. 

(e) The least upper bound of P[pn(En)c6’ | 6’’] converges to a withn > ~. 

(f) For any sequence of setfunctions {p,(E,)} which satisfies (d) and (e), the least 
upper bound of 


Plon(En)c6’ | 0”) — Plon(En)c6’ | 6”) 
with respect to 6’ and 6” converges to 0 withn > ~., 
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The validity of the above propositions follows easily from the identity 
Plon(E,)c6’ | 6] = 1 — P[W,(6’) | 0”). 
3. Assumptions on the probability density function. For any function 


¥(x) denote by Eef(x) the expected value of ¥(x) under the assumption that 6 
is the true value of the parameter, i.e. 


Ewa) = [ vaifle, 0) de. 


For any 2, for any positive 5, and for any real value 6’ denote by > 6’, 5) the 
greatest lower bound, and by ¢,(z, 6’, 5) the least upper bound ofS — pe 108 f(z, 6) 
in the interval 6 —i56 <6< @ + 6. 
Throughout this paper the following assumptions on f(z, 6) will be made: 
AssumPTIon I: The expectation Eg, = 5 18 f(x, 0’) ts a continuous function of 


6’ and 6, and for any pair of saemnae {0} and {6,} (n = 1, 2,---, ad inf.) 
for which 


lim Eg, ., log f(x, 6.) = 0 


lim (6, — 62) = 


T= 


Furthermore 
a 2 
Ew EB log f(a, 6 | 


2 
is a bounded function of 6’ and 6’’, and Eg |3 log f(z, | = d(@) has a positive 


lower bound. 

Assumption II: There exists a positive value ky such that the expectations 
Eyox(z, 0”, 6) and Ey-y2(zx, 6’, 5) are uniformly continuous functions of 6’, 6’ and 
5 where 6 takes only values for which |5| < ko. Furthermore it is assumed that 
Ee [oi(a, 0’, 6) (¢ = 1, 2) are bounded functions of 6’, 6’ and.d (| 5| < ko). 

AssumPTIon III: The relations 


C Slay Od [ SMe, 6) dx = 0 


hold. 
The above assumption means simply that we may differentiate with respect 
to @ under the integral sign. In fact 


+00 
f(x, 0) dx =1 
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identically in 6. Hence 


er er 
Gl, feOd= =] fle, 6) de =0. 
Differentiating under the integral sign, we obtain the relations in Assumption III. 
AssuMPTION IV: There exists a positive n such that 


a 2+n 
Ey E log f(z, | 
is a bounded function of 6. 


4. Some theorems. The assumptions on f(z, 6) made in this paper become 
identical with the assumptions I-IV formulated in a previous paper’ if a certain 
set w involved in those assumptions is put equal to the whole real axis 
(—«, +0). Hence we can make use of all results obtained in that paper 


putting w = (—«, +). Among others, the following statements have been 
proved there: 


> 1 @ 
(A) Denote >> aie 39 18 S(ta, 9) by yn(0, Z,) and let R,(6) be the region 
a=l 


defined by the inequality y,(@, E,) > An(@) where A,(6) is chosen such 
that P[R,(6) | 6] = 1 — a. Then for any sequence of regions {Z,(6)} for 
which P[Z,(@) | 6] = 1 — a, the least upper bound of 


P[Z,.(8’) | 6”) — P[Rn(8’) | 6] 


in the set 6’ > 6 converges to 0 with n > o. 

(B) Let S,(6) be the region defined by the inequality y,(@, EZ.) < B,(@) where 
B,(6) is defined such that P[S,(6) | 6] = 1 — a. Then for any sequence of 
regions {Z,(6)} for which P[Z,(6) | 6] = 1 — a, the least upper bound of 


P[Z,(6’) | 6] — P[S,(6’) | 6”) 


in the set 6’ < 6 converges to 0 with n > o. 

Denote by 7',(6) the region defined by | y,(@, EZ.) | > C.(6) where C,(6) 

is chosen such that 

(a) P(T.(6) | 6] = 1— a. 

Then T,,(@) satisfies also the following two conditions: 

(b) The greatest lower bound of P[T,(6’)6’’] converges to 1 — a with 
n— @®, 

(c) For any sequence of regions {Z,(6)} which satisfies (a) and (b), the 
least upper bound of 


P(Z,(8’) | 6”) — P[T.(6’) | 6”) 


converges to 0 with n > o~. 


(C 


a 


4A. Watp, ‘‘Some examples of asymptotically most powerful tests,’’ Annals of Math. 
Stats Vol. 12 (1941), pp. 396-408. 
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On account of Propositions I and II we easily get the following theorems: 
THEOREM I: Denote by £,(E,) the set of all values of 6 for which y,(0, En) < 

A,(6) and A,(6) is defined such that P{y,(@, En) > An(0)| 0] = 1— a. Then 

£,(EZ,,) satisfies the following two conditions: 

(a) Pl§.(EZ,)c0 | 6] = a for all values of 8. 

(b) For any sequence of setfunctions {t',(E,)} which satisfies the condition (a), 

the least upper bound of 


Plén(En)c6’ | 6] — Pl. (E,)c6’ | 0] 
in the set 6’ > 6 converges to0 withn > ~. 

THEOREM II: Denote by £,(E,) the set of all values of @ for which y,(0, E,) = 
B,(0) and B,(@) ts defined such that Ply,(6, E,) < B,(6)| 6] = 1— a. Then 
tn(E,) satisfies the following two conditions: 

(a) Plt,(E,)cé | 6] = a for all values of 6. 
(b) For any sequence of setfunctions {t,(En)} which satisfies the condition (a), 
the least upper bound of 
Pltn(En)c@’ | 6] — Pltn(En)c6’ | 0”) 
in the set 0’ < 0 converges to 0 withn > ~. 

THeEoreEM III: Denote by p»(E,) the set of all values of 0 for which | y,(0, E,) | < 
C,(@) and C,(8) is chosen such that P{| yn(0, E.) | > Cn(0)| 6] = 1 — a. Then 
pr(E,) satisfies the following three conditions: 

(a) Plon(E,)cé | 6] = a for all values of 0. 
(b) The least upper bound of P{p,(E,)cé’ | 6’’] converges to a withn > ~. 
(c) For any sequence of setfunctions {p,(En)} which satisfies the conditions (a) 
and (b), the least upper bound of 
Plon(E,)c6’ | 6’) — Plp.(E,)c6’ | 6] 
converges to zero withn — ~. 

Now we shall investigate the question whether the sets ¢,(Z,), ¢,(E,) and 
p.(E,) are intervals. For this purpose we will prove some propositions. 

Proposition III: Let ¢ and D be two positive numbers such thate < D. Denote 
by Q,(0, «, D) the region which consists of all points E, for which 


y(0+¢,E,) <—n', and = y,(@—€',E,) >n' 
for all values ¢’ in the interval [e, D]. Then we have 
(1) lim P[Q,(6, ¢, D)|@] = 1 


uniformly in @. 

Proor: Let «:, €2, °°: , €- be a sequence of points in the interval [e, D] such 
thatep -e = @ —a=-'+ = €— €& 1, = D — «€, = ky (Say), wherer is chosen 
sufficiently large such that Assumption II holds for |5| < ko. Denote by 
R,(6, ¢;) the region in which 


(2) yn(0 + 6, En) < —n'. 
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We will show that 
(3) lim P[R,(6, «) | 6] = 1 


uniformly in 6. 
From Assumption I it follows that the greatest lower bound of 


0 , 
Ey =, log f(z, 6 + é’) 


with regard to ¢’ in the interval [e, D] is positive. Let this greatest lower bound 
be A > 0. Since on account of Assumption I E, 5, log f(z, @ + e’) is a continu- 


ous function of ¢’, it does not change sign in the interval e < e’ < D. Since 


2 2 
this is true for arbitrarily small ¢ and since Ey, E log se.) | = —-E, = 


i log 


f(x, 6) has a positive lower bound (Assumption I), it follows easily on account of 
Assumption II that 


Ey © log f(z, 6+’) <0. 


Hence 
(4) Ey 5 log f(z, + ¢) < -A <0 for exe’ < D, 


and therefore 
(5) Ev.(@ + €, Ex) < —AV/n fore < e < D. 


From Assumption II it follows that the variance of y,(@ + ¢’, E,) is a bounded 
function of @ and e’. Hence 


(6) lim Plyn(@ + 6, En) < —3AVn|6] = 1 


uniformly in 6. The equation (3) is a consequence of (6). 
Denote by S,(@, ¢;) the region in which 


> Delta, @+6,b)|<C @=1,2) 


where C is greater than the least upper bound of | Eyg,(z, 6’, ko) | with respect 
to 6and 6’. Then we have on account of Assumption IT: 


(7) lim PIS, |= 1 = 1,2,---,9) 


uniformly in @. In the region S,(6, ¢;) we obviously have 


(8) Yn(O + €:, En) < ya(O + €&, En) + 2kov/nC 
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for all values ¢; in the interval [e; — ko, €: + ko]. By choosing r sufficiently 
large we can always achieve that 


2hC < 


| pe 


Denote by T’,,(8, €:) the region in which 
9) wl +6, Bs) <—S Vn for a— bb SiS ath. 


From (6), (7) and (8) we get 
(10) lim P[T,(6, «) |] = 1 


uniformly in 6. Let Q/,(6, «, D) be the common part of the r regions 
T,.(0, &), °°: , Tn(9, €), i.e. Q,(8, €, D) is the set of all points EZ, for which 





yal0 + €, B,) <—4 Vn 


for all ¢’ in the interval [e, D]. 
on n, we get from (10) 


(11) lim P[Q;.(6, «, D)|6] = 1 





Since r is a fixed positive integer not depending 


uniformly in 6. 
In the same way we can prove that 


(12) lim P[Q‘.(6, «, D) | 6] = 1 


uniformly in 6, where Q’,(6, ¢, D) denotes the region in which 


yn(@ — é, E,) => : a/n forall ¢’ in [e, D]. 


Proposition III follows from (11) and (12). 
Proposition IV: Denote by V,(0, €) the region in which 


a , $ 
30 yn(0 ’ E,) ~ -@ 


for all values 6’ in the interval [@ — e,9+ «]. There exists a positive ¢ such that 
lim P[V,(@, «) | 6] = 1 
untformly in 6. 


Proor: Since the least upper bound of Eeg.(z, 6, 0) is <0, we get from 
Assumption II that the least upper bound of Eey.(z, 6, €) is <0 for sufficiently 
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small e > 0. Denote the least upper bound of Egge(z, 6, €) by —B and let the 
region in which 


~ © ealte, 3,6) < 3B 


be denoted by W,,(6, €). From Assumption II it follows that 
lim P[W,(6, «—)|@] = 1 


uniformly in 6. Since for almost all n W,,(@, €) is a subset of V,(6, €), Proposi- 
tion IV is proved. 

Proposition V: Let A,(0), Bn(6), Cn(0) be the functions as defined in Theorems 
I-III. There exists a finite value G such that 


| A,(0) | < G, |B,(0)|<G and |C,(0)| <G 

for all 6 and all n. 

Proposition V follows easily from the fact that the variance of y,(6, E,) is a 
bounded function of ~ and @. 

Let D be an arbitrary positive number and denote by W,(6, D) the region con- 
sisting of all points Z, for which the following conditions are fulfilled: 

(a) The equation y,(6’, E,) = A,(6’) has exactly one root in 6’ which lies in the 
interval [@ — D, 6+ D). 

(b) The equation y,(@’, E,) = B,(6’) has exactly one root in @’ which lies in 
the interval [@ — D, 6 + D}. 

(c) The equation y,(6’, E,) = C,(6’) has exactly one root in @ which lies in 
the interval [(@ — D, 6+ D]. 

(d) The equation y,(6’, E,) = —C,(6’) has exactly one root in 6’ which lies in 
the interval [6 — D, 6 + D]. 

(e) The common part of [@ — D, 6 + D] and ,(E,) is the interval [6;,(E,), D] 
where 6,,(E,) denotes the root of the equation in (a). 

(f) The common part of ¢,(E,) and [@ — D, 6 +.D] is the interval [—D, 6,(E,)] 
where 6,,(E,) denotes the root of the equation in (b). 

(g) The common part of p,(Z,) and [@ — D, @ + Dj is the interval 
[0.(E,), 6.(E,)] where 0,(E,) denotes the root of the equation in (c) and 
6,(E,) denotes the root of the equation in (d). 

From Propositions III-V follows easily the following 
Proposition VI: For any positive value D 
lim P[W,(6, D) | 6] = 1, 

uniformly in 0, provided that the functions A,(6), B,(0) and C,,(6) are continuous 

and of bounded variation in any finite interval. 


We will show that Proposition VI remains valid for D = + ©, if we make the 
following 
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Assumetion V: Denote by ¥(z, 0, D) the least upper bound of x log f(z, 6’) 
with respect to 0 where & > 6+ D. Denote furthermore by y*(zx, 0, D) the greatest 
lower bound of 5 log f(z,6’) with respect to & where 6’ < @— D. There exists a 


positive D such that the least upper bound of Ew)(x, @, D) with respect to 0 is negative, 
the greatest lower bound of Ew)*(x, 0, D) with respect to 0 is positive, and the variances 
of v(x, 0, D) and y*(z, 6, D) are bounded functions of 6. (The variances are 
calculated under the assumption that 0 is the true value of the parameter.) 

It follows easily from Assumption V that 


lim P| F-Lvle 6, D) en ‘0 


= lim P| FLV es D) > nto =1 


uniformly in @. 
Since 


FL Vtu, 0, D) > yal", Ba) for > 0+ D 


and 


a > ¥*(xa, 9, D) < yn(O’, Ex) for 6 <0—D, 


Proposition VI remains valid if we substitute + for D. 
Hence we obtain the following 
Corotuary: If the assumptions I-V are fulfilled and if A,(@), B,(0) and C,(6) 
are continuous and of bounded variation in any finite interval, then 
(a) The root 6,,(E,) of the equation y,(0, E,) = An(6) in 0 is an asymptotically 
best lower estimate of 0. 
(b) The root 6,,(E,) of the equation y,(0, E,) = B,(6) in 6 is an asymptotically 
best upper estimate of 0. 
(c) The interval (8,(E,), 6.(En) 1s an asymptotically shortest unbiased confidence 
interval of 6, where 0,(E,) denotes the root of the equation y,(0, E,) = +C,(@), 
and 6,(E,) denotes the root of the equation y,(0, E,) = —C,(8). 


5. Some Remarks. 1. I should like to make a few remarks about the relation- 
ship of these results to those obtained by S. S. Wilks.® The definition of a 
shortest confidence interval underlying Wilks’ investigations is somewhat differ- 
ent from that of Neyman’s which has been used in this paper. According to 
Wilks, a confidence interval 5(E,) is called shortest in the average if the expected 


5S. S. Witxs, ‘‘Shortest average confidence intervals from large samples,’’? Annals of 
Math. Stat., Vol. 9 (1938), pp. 166-175. 
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value of the length of 6(Z,) isa minimum. The main result obtained by Wilks 
can be formulated as follows: The confidence interval [@,(E,), 6.(Z,)] given in 
our Corollary is asymptotically shortest in the average compared with all confi- 
dence intervals computed on the basis of functions belonging to a certain 
class C. In the present paper no restriction to a certain class of functions has 
been made. 

2. If the parameter space 2 is not the whole real axis, but an open subset of 
it, and if the assumptions I-V are fulfilled when 6 can take only values in Q, 
the previously proved Corollary remains valid. If 2 is a bounded set, Assump- . 
tion V is a consequence of Assumptions I-IV. 











GROUPING METHODS 
By P. S. Dwyer 


University of Michigan 


1. Introduction. The conventional formulas for moment adjustments known 
as Sheppard's corrections are not too satisfactory for practical use. As Carver 
has pointed out [1] Sheppard’s corrections are merely systematic adjustments 
which eliminate the bias introduced by grouping. The values of the moments 
after Sheppard’s corrections have been applied may be looked upon as unbiased 
grouping estimates of the true moments while the uncorrected values constitute 
biased estimates. 

In practice one obtains his moments from a single grouping. The application 
of Sheppard’s adjustments in such a case does not necessarily result in the un- 
biased estimate being closer to the true moment than is the biased estimate and, 
in an appreciable percentage of cases, the unbiased estimate is further from the 
true moment than is the biased estimate. One does not know when he applies 
Sheppard’s adjustments to the results of a single grouping whether or not he is 
making a correction in the right direction. 

This situation is not too satisfactory and yet practical necessity demands 
some method of grouping. The improvement of modern calculating machines 
tends to push grouping techniques further into the background since, in many 
cases, the machines permit the determination of the actual values of the moments 
without grouping in a reasonable amount of time. But even here it is possible 
to use grouping methods and to get a good estimate of the true value in a frac- 
tion of the time. It is the purpose of this paper to present some new grouping 
methods which are useful in obtaining much better unbiased estimates from a 
single grouping than can be obtained with the use of Sheppard’s corrections. 
These methods demand additional work but this additional work is justified by 
the additional precision resulting when such precision is desired. 

The spirit of the new approach, which in one sense is a generalization of the 
earlier approach, can be expressed very simply though the details of the develop- 
ment and the calculational methods demand amplification. If we let « = the 
true value and 2’ the grouped value (the value of the class mark of the group in 
which x is), and the error, « = the difference between the true value and the 
grouped value, then ° 


(1) e=2z-7', r=2’+e, and v’=2r-e«. 


ss 


In the classical theory we use Sx” as the biased grouping estimate of Z2*. In 
the new methods we use Yxx’*' as the biased estimate or, if we desire more pre- 
cision, 2°’ * as the biased estimate. It is then possible to correct Zxz’*™ for 
grouping bias and to correct Sa’x” for grouping bias just as we now correct 
Sx" for grouping bias. It is also possible to use the values of =z” and Srz’** 
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in obtaining a better unbiased estimate or to use the values of Ez”, Ezzx'*", 
and >2’x"*” in obtaining a still better unbiased estimate of Dz". 


2. Illustration. The relative merits of the conventional method and the 
proposed methods can be shown effectively by means of an illustration. For 
this purpose I have selected the problem used previously [1, 154] in showing the 
variations in grouped results. The power sums rather than the moments are 
used and the origin is taken at a point near the mean so that the relative varia- 
tions are as large as possible. If the values of the power sums were “padded” 
by measurement about zero, the relative variations would not appear as large. 
However, a problem which shows considerable variation, and in the problem 
under consideration the nine =x” unbiased grouping estimates of 22° resulting 
from the nine groupings do not even have the same sign, is an appropriate one 
with which to demonstrate the improvements introduced by the new methods. 

The problem consists of 244 discrete variates which range in value from 64 
to 155. Carver took a class interval of nine and formed the nine frequency 
distributions which result when class intervals of nine are chosen in all possible 
ways. He computed the values of D2’, Dx”, Dx", Dx", for each of the nine 
distributions, corrected each for bias with the use of the Sheppard adjustments, 
and showed that the averages of the nine corrected estimates are respectively 
the values of =z, 22’, E2°, U2". 

In Table I are presented the values of the biased and unbiased grouping 
estimates of Zax, Da’, Zx*, Zax‘, which result from the use of (1) Dx"; (2) Saar’; 
(3) Sa?x'**; (4) Da", Dra’; and (5) Tx", Vrx’*", Vax’. The results are 
presented here for comparison only; the details of the computation are explained 
later. Rows of biased estimates are indicated by B while the rows of unbiased 
estimates are indicated by U. Parentheses are used to indicate entries which, 
while appearing in rows of biased estimates, are actually unbiased. The exact 
values of 2x", when they appear, are indicated by underscoring. The Roman 
numerals indicate the different frequency distributions while the grouping 
methods are indicated by the values of (1), (2), (3), (4), (5) above. The true 
values are 2x = —129, D2” = 77,591, Sa*° = —52,005, Dx* = 69,239,951 where the 
values of x used are the values of the original variates decreased by 105. 

The information contained in Table I deserves more than cursory examination. 
Study shows that the estimates resulting from method 2 are much closer to the 
true values than are the estimates resulting from method 1,-etc. Table II is 
presented below in order to facilitate the comparison of the relative amounts of 
grouping error involved in the different methods. The standard deviation of the 
grouping error of the conventional method, method 1, is used as a norm and the 
standard deviations of the grouping errors for the new methods are compared 
with this norm. 

The decline in the size of the error revealed in Table II indicates a decided 
decrease in grouping errors. Grouping method 2 enables one to compute the 
mean exactly and this is always possible when method 2 is applied to discrete 
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TABLE I 






Biased and Unbiased Grouping Estimates by Different Methods 


Grouping 
Method 


1-B 
2-B 
3-B 
1-U 
2-U 
3-U 
4-U 


2-B 
3-B 
1-U 
2-U 
3-U 
4-U 
5-U 


2-B 
3-B 
1-U 
2-U 
3-U 
4-U 


2-B 
3-B 
1-U 
2-U 
3-U 
4-U 
5-U 


77,591 


(77,591) 
75,5224 
76 ,593 


77,591 
77 ,6634 
77,591 








78 , 466 
(77,181) 


(77,591) 


76,839} 
77,181 
77,591 
77,5223 
77,591 





77,769 
(76,797) 
(77,591) 
76,1424 
76,797 
77,591 
77,4513 
77,591 








79,747 
(77 ,790) 
(77 ,591) 
78,1204 
77,790 
77,591 
77 ,459% 
77,501 








81,934 
(78,891) 
(77,591) 
80,3074 
78,891 

77,591 

77,4744 
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—54,602 
—52,977 
(—52,465) 
—50, 242 
—52,117 
—52,465 
—52,307 
—53,066 


2,889 
—17,037 
(—35,097) 
5,109 
—16,177 
—35,097 
—59,469 
—51,291 


—23,311 
—34,464 
(—44, 108) 
—20,531 
—33,604 
— 44,108 
—59, 350 
—52, 243 


19, 666 
—4,521 
(—28, 387) 
21,746 
—3,661 
— 28,387 
—55,475 
—51,932 








zx 
69, 239,951 


69 , 063 , 265 
68 , 267 ,577 
68 ,033 ,657 
66,023,177 
66,735,717 
67 ,516 , 383} 
69,001,817 
69, 190 , 537 


74,519,962 
72,296,367 
70,685,801 
71,427,194 
70,752,747 
70, 168, 5273 
68,770, 406 
69, 246, 172 





71,465,409 
70 ,053 ,093 
69, 214,545 
68 , 400 , 521 
68 , 517, 153 
68 , 697 , 2713 
68 , 945, 609 
69 , 201 , 169 





74,171,443 
72,095 , 664 
70,592,774 
71,027 , 435 
70 , 539 , 864 
70 ,075 , 500 
69 ,037 ,511 
69 , 248 ,077 








76,143,874 
73,590,207 
71,610,053 
72,912,386 
72,012,387 
71,092,7793 
69, 142,430 
69,312,700 
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zx 
69, 239,951 


72,467 , 541 
70,940,124 
69 , 902,910 
69 , 307 ,613 
69 , 379 , 524 
69 , 385 , 636% 
69 , 536 , 657 
69 , 241,507 


71,851,930 
2-B 70,515,354 
3-B 69, 647,462 
1-U 68, 685,722 
2-U —129 68,951,994 
3-U = 69, 130, 1884 
69 , 689, 930 


4-U —129 
5-U —129 69 , 260, 146 


1-B (—89) 
2-B (—129) 
3-B — 
1-U —89 
2-U =129 
3-U _ 
4-U —129 


5-U “  —129 


68,426,497 
67 959,816 
67,944,476 
65,330, 249 
66,412,776 
67 ,427 , 2023 
69,711,437 
69, 210, 235 


—82,416 
(—65,788) 
—99,577 
—81,556 
—65,788 
—47,114 
—51,473 


1-B 
2-B 
3-B 
1-U 


(—180) 
(—129) 


—180 


78 ,894 
(77,517) 
(77,591) 

77 , 2673 


— 180,792 
— 134,865 
—177,192 


73,155,150 
71 ,407 ,737 
70, 183,341 
70 , 045 , 262 


2-U —129 77,517 


3-U - 77,591 
4-U —129 77 , 7664 


5-U —129 77,591 


—134,005 
—92,169 
— 45,591 
— 52,704 


69, 857,397 
69, 666, 6664 
69, 323, 762 
69, 246,016 


data. There is also a corresponding decrease in the errors of the higher powers 
to roughly one-half, two-thirds, three-fourths. Greater precision in the case 
of the higher power sums can be obtained with the use of the other methods, 
though these methods demand more calculation. 


There is one more question which should be discussed before the general 
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theory is presented, and that deals with the method of computation of the 
quantities =zz’*", Da°x’**, in methods 2 and 3. Computational techniques are 
discussed in a later section of the paper, but enough should be given now to 
make the meaning of Drx’ and 2272" clear. In getting =x”, we recall, 
we need only the values of the class mark, x’ and the frequency associated with 
each, f,. To get the values of =zx’** we need in addition to x’ the sum of 
the x values which are grouped together in the class having the same class mark, 
x’. We denote this value by x, and we use this instead of the f. of the usual 
method. In the case of method 3 we record 2x2, where 22, is the sum of the squares 
of all x values having the same grouped value z’. 

Let us examine the first grouping in Table I. The original 244 variates were 
recorded by Carver [1, 154] and he gave the values of f, for each grouping. 
It is necessary for us to return to these original variates, but instead of counting 
the variates in a given group, we add them and we add their squares. 

In obtaining the values for the first grouping in Table I the variates were 
transformed with the use of x = v — 105. The variates then ranged from 


TABLE II 


Standard Deviations of the Grouping Errors of the Different Methods Expressed as Percentages 
of the Standard Deviations of the Usual Method 





Method =x 











| =x? | =x | zx 
1 100 | 100 | 100 | 100 
2 0 48.9 | 65.5 | 74.3 
3 - 0 32.1 49.1 
4 0 8.8 | 7.3 13.8 
5 0 | 0 | 9 1.5 
—41 to 50 and the frequency distribution was made with mid values 2’ = —37, 


—28, —19, —10, etc. The values of f., , x, , and x2 were then computed and 
recorded in the columns, 2, 3, 4, of Table III. The next three columns are 
computational columns useful in obtaining the biased estimates recorded at the 
bottom of Table III and also in Table I with the use of i Sse 


z 
2. re = - ten", ; gr? = z. real, 
= 3’ 


3. General formulas for corrections for grouping bias. We are next led to 
the question of correcting these estimates of 2x’ for the bias introduced by 
grouping. Before indicating the numerical work, we derive general formulas 
for correction for grouping bias. 

We assume that the variates are recorded in units of h which means that, in 
the case of discrete series, the smallest possible difference between any two 
unequal variates is equal to h. In case the distribution is continuous, the 
recorded values constitute a discrete series recorded in unitsof h. Thus heights 
may be recorded to the nearest inch, in which case h is one inch, or to the nearest 
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one hundredth of an inch, in which case h is 1/100 inch, etc. We assume 
further that all possible groupings of k different values are made. Thus if the 
smallest variate is a, then the values of x = a,a +h,a-+ 2h,---,a+k + th, 

-,a+k — 1 hare thrown ina group with class mark a + 4(k —1)h. The 
k; possible sets of groupings of k are made in this way. 

We examine the error involved when a specific variate x is replaced by the 
class mark zx’ in each of these groupings. The values of the lower open limit, L, 
the upper open limit, U, the class mark, x’, and the error e = x — 2’ are in- 
dicated in Table IV. The k different groupings indicated by the different rows 
show zx at the lower limit, x one step above the lower limit, x two steps above 


TABLE III 


Values of x’, fe’, xz: and x*,’ for the First Grouping with Computation of Biased Estimates 
of =x" 


148,877 | 7,890,481 
85,184 | 3,748,096 

42,875 | 1,500,625 

17,576 456 , 976 

4,913 83,521 

4,096 

-1 1 

100 | 1000 10,000 

! | 361 —6859 | 130,321 
784 «| —21,952 614 ,656 

1369 | —50,653 | 1,874,161 


—134,191 | 68,267,577 | 
69,063,265 | | 


the lower limit, etc. It is at once apparent that the errors in replacing 2’ for z in 
the k different ways constitute the deviations from the mean of the rectangular 
distribution h, 2h, 3h, --- , kh. We indicate the moments about the mean of 
this rectangular distribution by Ri, R:, R;, R, and we use the notation E(¢') 
as the sum of the ith powers of the k different e’s divided by k. It follows that 
E(e') = R,. Now the values of R; are 0 when t¢ is odd and are well known when 
t is even [2,325]. The ones in which we are especially interested are 


_kK-1..2 _ (kK — 18k — 7) 14 
me a & and eee A 
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If an adjustment of scale is made so that the differences between successive 
class marks are unity, as is customary, the value of h is 1/k. The values of 
R, and R, are then 

_1-1/k _ (1 — 1/K)(3 — 7/F) 
” EE e—— 


As the number of groupings increases the value of 1/k” > 0 and the appropriate 
values of the moments of the continuous rectangular distribution result. Thus 
1 7 


—_ a 1 _ = 
Re = 55: fu = 5» and 6R2 Ru = 990° 







































TABLE IV 
Open Limits, Class Marks and Errors for the Different Groupings 





x’ = 4(L + U) e=xz— x’ 


—4(k—-1)h 


i | | a | | 


—3(k—3)h 


ec | ef SS SSS sSSSS 


—4(k—5)h 


— ee ee 


| | a | | 


—4[k—(2i—1)]h 


SS 


NESS 


4[k—3]h 


Oe 


$[k—1)h 


If now we let F, be any real function of x defined for the values x = a, a + k, 
a + 2k, --- , we have at once the useful lemma 


(4) E(zz* é F,] = 22°F.E[e'] = R:rz’F. . 


This results from the fact that the values of x, and of all functions of z, are 
unchanged by the groupings even though the values of x’ and e¢ vary. 

The = in (4) indicates a summation with respect to the variates while the 
summation with respect to the different errors is taken care of in the E notation. 
The limits of the = in (4) are purposely left indefinite so that either a serial 
or @ frequency notation can be used. Thus if a serial notation is used, the 
limits are from 1 to N, the values of z are the variates z; and F, becomes F,, . 
In this case F,, may be set equal to unity to give the corrections of method 1, 
may be set equal to 2; to give the corrections of method 2, or may be set equal 
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to x; to give the corrections of method 3. In case the notation of the frequency 
distribution is preferred, the limits of the summation are the smallest variate 
and the largest variate, the values of x are the values of the different variates 
which occur. In this case we may have F, = f, , the frequency function, 
F, = af» = 22,01 F, = 2 fe = x2. 

The continued application of (4) to the terms in the expansion of E[=z’F,] 
results in 


E[> x" F.] = E[> (« — ©'F.] = E| > (—1)' () oF, | 
=o) dere = Ev) aden. 


The fact that R; = 0 when ¢ is odd may be used in writing out the expansion. 
It is possible to work out a more general theory where the class mark is some 
other value (say the smallest variate) rather than the mid-value. In such a 
case formula (5) would apply, but the values of R, would be the values of the 
moments of a rectangular distribution rather than the central moments. The 
above formula is sufficiently general for the purposes of this paper. 

Specific values of (5) when s = 0, 1, 2, 3, 4 are 


E(=F.] = =F, 
E|Z22’F,|] = 22xF, 
E(3z"F,] = 22°F. + R-=F; 
E(=z"F,) = 22°F, + 3R2=2F, 
E(=z"F.] = =2‘F, + 6R.=a°F, + Ri rF-.. 


(5) 


These equations can be solved for 2F,, =xF,, etc., in terms of the expected 
values. If we use the inverse operator and write E'[B] = A instead of E[A] 
= B we have 


E™[=F,] = 
E™[2F,] = 
E"(32°F,] = =2"F, — R2=F, 
E“|[2°F.] = 22°F, — 3R.=2'F, 
E"[>a‘F.] = =2"F, — 6R.2x"F, + (6R} — R,)=F; and in general, 


E"[>> z'F.] = > x"F. — > (—1)'()R.E" [> x* ‘F.]. 


These values E~'[=z'F,] are unbiased estimates of =2'F, since 


E(E™[22'F.]] = =x'F, 
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The corrections for method 1, the customary corrections, are obtained if a 
serial notation is used with F, = 1. The corrections for method 2 are obtained 
if a serial notation is used with F, = x. The corrections for method 3 are 
obtair d with F, = xz. Thus we have 


E"(N] = N 
E™ [2] = 32’ 
(8) E™ [227] = =z” — RN 


E™ [22°] = 22" — 3R, 2’ 

E™[=2‘] = 22 — 6R.=x” + (6R: — RN; 
E“ [=z] = 22 

E™ [22°] = Zaz’ 


(9) E™ [32°] = 2axr” — R.rx 

E™[22‘] = 22x" — 3R.22xx’ etc., 
and 

E™ [27] = 22° 
(10) E™ [32°] = 22°’ 


E™[2‘] = Sa’x” — R.=z’. 
These formulas are the ones used in obtaining the unbiased estimates in methods 
2 
1, 2, 3 from the biased estimates in Table I. In this case R2 = oS ee 


2-3 
2 $ 
6R2 — R, = = - 2 0 ~ Sak wed the whe fellow ty 


direct substitution in (8), (9), (10) above. 





















4. Compound grouping formulas. So far nothing has been said about the 
calculation of the results by methods 4 and 5. These methods might be called 
compound grouping methods, since they utilize the biased results of more than 
one grouping method. The values of =x’ and Sxrz’*” are needed for method 4 
and the values of 2x", Sax’, Dx*x'*~* for method 5. The formulas for method 
4 are first presented. The argument is given in some detail for the value of 
E™[22’]. Now 


Da? = Yr! + 6)? = De” + Wa'e + TE 
= Ye? + 23a'(4 — 2’) + Te 


= — Dz" + 2Qdrr' + Te. 
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If the values of ¢ are known, we would have the exact value of 22’ since we know 
>a” and Zxx’. However, we do not know these values of ¢ from a single group- . 
ing, so we derive a formula giving unbiased estimates of Zz’. We have at once 


E[z2*] = 2a? = E[—Zx” + 2222’) + NR 
and since NR, = E[NR,] we have 
Sa? = E[— Dx” + 222x' + NR} and 
E“[27] = —Za2” + 222r' +-NR2. 


There is a relatively small error in this estimate since the only error involved 
is the difference between NR, and the actual sum of the squares of the e’s. 
This formula is the basis of the values of E™'[Zz’] recorded in Table I under 
method 4. For example, in grouping I, the estimate is —77149 + 2(76593) + 
244($8) = 77663% and this differs by only 723 from the exact value. 

In a corresponding manner, we may prove 


E“[22°] = —222" + 3222" + 3R.rz 
E22] = —332" + 4222" + 6R.E™[22"] + 3NRy. 


Different values of E~'[2z*] can be used. Im the calculations of Table I the 
values E[2"] = Zaz’ from (9) were used, but the values E[E27] = — 2a” 
+ 222’ + NR, could be used to give somewhat better results. 

It can be shown also that 


(11) 


(12) 


(13) E“[sa°] = —422" + 522r" + 10R,E™[22*] + 15R,zz; 
E“ [32°] = —52za’® + 622r" + 15R.E [Zr] + 45R,.E [227] + 5RN, 


and, after some argument that 


ES 2] =-@-1) Dae" +s yD 2""s 
14 te] : 
(14) >> (, (2t — 1)RxE [SD 24] 


where [4s] indicates the integer $s or $(s — 1). 
It is possible to obtain better unbiased estimates if we use in addition the 


values of 2a*x’"*. In this case the values of Zz and =z’ are known exactly, 


and after expansion of 2(x’ + «)’, replacement of e by x — x’ and of ¢ by (x — 2’)’, 
and further reduction, we get 


E“(d2*] = 22" — 32ax" + 3E2°2’, 

E22] = 322" — 8222" + 622°” — 3NR,, 

E” [2] = 622" — 15222" + 1022°2”" — 15R,2z, 

E™[d2°] = 102° — 24222" + 1522°2" — 45R,22* — 10NRe, 


(15) 
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- and in general, 


E* [> z'] = #(s — 1)(s — 2) > 2” — 3(s — 2) Dar” 
(16) [ts] om 
+ 38(s — 1) a et? — 2. (3,)(7 ') Ry, E™* (dX a *4), 


t=_2 2 2 


Compound formulas involving additional quantities such as 2a°x”*, Ea‘x’**, 
etc., can be worked out by the methods outlined above. 
5. Computational methods. It has been shown in sections 3 and 4 how the 
unbiased estimates can be obtained from the biased estimates. It is the purpose 
of the present section to show how these biased estimates can be computed 
efficiently. One method of calculation was shown in Table III. The values 
of fe, 2, 22 were computed and recorded, and the resulting power sums 
obtained. This is the most direct means of computation and if the number 
of groups is small and if a modern computing machine equipped with auto- 
matic positive and negative multiplication is available, it may be the preferred 
method. It should be noted that the values of 22, in Table III are obtained 
most easily with the use of a machine which permits the calculation of the 
square with a single key punching operation. 

It is customary to use the devices of subtraction of a constant (either a central 
class mark or the smallest class mark) and division by a constant (size of the 
class interval) to simplify the computational work. Thus in Table III we 

/ 
—— and compute the_values of 
Yd"F, . If Fz is the frequency function, we have the usual formulas, but if 
F, is x, or xz , then the results are terms of the type 22d’ or Za'd""*. It is 
x — (-—37) 

9 





could use the transformation d’ = 


possible to reduce these to equivalent variables by the use of d = 


so that the values of 2d”, Sdd’*", =d’d’* result. We then correct for bias 
with the use of the formulas of sections 3 and 4 where the power sums of the 
rectangular distribution are computed with h = 1/k. 

Another method which in many cases is preferable to that just described is 
the method of cumulative totals. The values of f,, , z2 , and 2, , are cumulated 
successively for the different values of x’ and the values of the biased grouping 
estimates are obtained immediately from the entries in the last few rows. The 
cumulations of Table III are shown in Table V. The entries in the column 
of the highest cumulations of f, , zz, x2, with the exception of those at the 
bottom of the column, need not be recorded. 

It is possible to provide multipliers for these entries by an adaptation of 
a method given in an earlier paper [3]. A table of multipliers has a top marginal 
row composed of a, a + k, a + 2k, etc., and a left marginal column composed of 
k — a,2k —a,ete. The first row in the table is composed of 1, k — a, (k — a)’, 
(k — a)’, etc., and the first column of 1, a, a’,.a’, etc. Each entry in the table 
is found by adding the product of the entry above it and the columnar heading 
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to the product of the entry at the left and the row heading. The multipliers 
for a = —37 are shown in Table VI. 

The diagonal terms are the multipliers of the values of a given cumulation. 
Thus the multipliers of the bottom entries of the columns of each of the three 
sets of cumulations of Table V are successively 1; —37, 46; 1369, —3323, 
2116; ete. 

This method is ideally adapted to the use of Hollerith cards. The information 
is punched on the cards to the number of places desired. The computational 
grouping is then accomplished by sorting. As an illustration we take the 








TABLE VI 
Multipliers when a = —37 and k = 9 
—3? | —28 | -19 | —10 wf 
46 | 1 | 46 | 2,116 | 97,336 4,477,456 
55 | —37 | —3,323 | — 222,969 — 13,236,655 
64 | 1,369 | 180,660 | 15,798,651 | 
| 


73 | — 50,653 —8,756, 149 
82 | 1,874,161 | 


TABLE VII 
Hollerith Illustration 




















x” = x! —4.5 | CU.) | C(xz,) | C(yz+) 
210 2 | 422 | 1,387 
200 | 6 | 1,242 | 4,153 
190 10 | 2,017 | 6,961 
180 21 | 4,035 | 14,733 
170 48 | 8,727 | 33,675 
160 110 | 18,923 | 76,565 
150 250 40, 466 | 173,491 
140 458 | 70,392 | 316,073 . 
130 719 105,434 | 492,774 
120 | 900 | 127 , 990 | 613,763 
110 | 980 137 , 203 666 ,088 
100 | 999 | 139, 199 678, 299 

90 | 999 | 139, 199 678 , 299 
80 | 1,000 } 139, 288 678, 896 








records of the weights of 1000 students as reported by Carver [4] when measured 
to the nearest pound. The value of =z is 139288 lbs. and that of =z” is 19,692,450 
(Ib.)’ and we wish to obtain approximations to these values by grouping. If we 
let the grouping intervals be 80-89, 90-99, etc., with class marks 2’ = 84.5, 
94.5, 104.5, etc., we would find by usual methods Sx’ = 139,520 lbs. and Sz” = 
19,760,430 (lb.)”. However, it is possible to wire in the three place number z, 
and to get from the same number of groups =x = 139,288 lbs. and Yzz’ = 
19,727,326 (Ib.)”. The unbiased values for method (1), (2), or (4) can be com- 
puted with the appropriate formulas of sections (3) and (4). 
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The Hollerith run is shown in Table VII where the first column indicates 
the smallest variate in the class rather than the class mark. The next columns 
show C(f,) and C(z,). Thefourth column C(y,-) is discussed in a later section. 

The values for method 3 and method 5 cannot be obtained so readily, since 
the quantities to be grouped are the 2° and these do not appear on the card. 
However, it is possible to use a multiplying punch or to use a table of squares 
in the form of prepunched cards to get these values of x” on the cards. It might 


be preferable, in some cases, to do this work and then to use a coarser grouping 
than would be used otherwise. 


6. Moments. The formulas (7), (8), (9), (10), give moment formulas if the 
w./P ws ,./P.,.9 


. + alk — Zz 
proper values of F, are assigned. We let », = —— and »,, = —,— and have, 


N N 
in case F, = 1/N in (7) the usual formulas 


= 
=n— kz 
v3 — 3k 
= vw — 6m. + (GR: — R,). 
If F, = xr/N we have 
E™ [ui] = M1 
E [us] — waa 
E™ [us] = vn — Row 
E™ [us] = v3 — 3Rovn . 
While if F, = 2°/N we have 
E™[us} = we 
(19) E™ [us] = ae 
E™ [us] = vee — Rowe. 


(18) 


Similar formulas can be written for methods (4) and (5). 
Previous to Carver’s article in 1936 it was assumed that central moments 
could be used in place of moments in formulas (17) without introducing bias, 


but this article demonstrated that estimates obtained in this way are slightly 
biased. Thus 


= E(v. — vj) = E(m) — E(v}) 
(20) = po + Re — po(v1) = we + Re — [fie(nn) + ui] 
= fio + Re — jio(v) so that 
E™ [jie] = 5. — Re + fia(r1) 


and so % — Rz is a biased estimate of je . 
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The general question of unbiased estimates of the central power sums and 
the central moments is one which has been studied for the conventional case by 
Pierce [3] and Craig [5]. The more general discussion resulting from the intro- 
duction of the new methods is one which may well be deferred to a later paper. 
It is interesting to note in passing that the estimate of the variance obtained by 
substituting central moments for moments in method (2) is not biased since 


E(on) = Ely — nym) = we — Mi = fic. 


It is to be noted that the formulas previously used give correct results 
when the adjustments are defined to make the power sums and the moments 
rather than the central power sums and the central moments unbiased with 
respect to grouping. A sensible method of procedure in such a case is to make 
the correction on the power sum as soon as it is computed. 


7. Product moments. Correlation. The introduction of additional variables 
opens up a variety of situations, since each of the variables may be grouped in 
different ways. Of these situations, one is immediately solved with the use 
of the formulas of section 3, and that is the case when one of the variables is 
not grouped. Let y be the ungrouped variable and let F, = y, be the sum of 
all the values of y having the same x grouped value, x’. This situation is 
frequently encountered when using Hollerith cards, as it is only necessary to 
wire in the whole variable y and take totals when the smallest value of z in the 
group is attained. Thus in Table VII, the value of C(y,-) can be obtained 
simultaneously with the value of C(f.) and C(z,-). Additional cumulation 
C(zz), C(wz-), etc., could be obtained at the same time. It follows from Table 
VII that 

Sy = 678,896 
E“3[xy] = 22x’'y = 94,929,322. 
The actual value of Lzy is 94,774,336. 

The general development of the theory of unbiased estimates of product 
moments is too extensive to be inserted here, but a brief outline might be in- 
dicated. We let the grouping errors be e = x — x’ andn = y — y’. Thenthe 
generalization of the lemma (4) is 


(21) > d ate yn’ F.G, = Rw Roa 22’ F: > y' Gy; 
z v z v 


where Ryo is the bth central moment of the rectangular distribution consisting 
of e’s and Rog is the dth central moment of the rectangular distribution con- 
sisting of n’s. This is applied in turn to the terms of the expansion of 


rx" y"F Gy . 
For example 


E{=2'y'F.G,) = Ex(x — ey — n)F.G, 
LryF G, = Ry XyF Gy, = Ru rzF.G, + RyRakF.G, ; 


(22) 








i 
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and if F, = 1 and G, = 1, we have 
(23) E{=2'y'] = Szy sothat E [Ezy] = =2’y’. 


If we use the customary device of correcting the moments for bias, rather 
than the central moments or the ratio which is the correlation coefficient, we 
have the usual formula for correction of the correlation coefficient in which 
the numerator term is not corrected for bias, but the values in the denominator 
are corrected. 

The use of method 2 gives =2’y’ and =zy’ as unbiased estimates of Lzy. 
It has been pointed out that these quantities 2x’y and rzy’ are readily obtained 
when the actual values of x and y are punched on Hollerith cards. Each is 
in general a better estimate of Zzry than is =2’y’ since one of the values in the 
product, in each case, involves no approximation. An average of these might 
be taken to obtain a better estimate of Zzy. If the values of =x’ and Ly’ are 
also available, it is preferable to use the formula 


(24) E“{r| = // Are Are where A,, = NExy — (Z2)(Zy). 
z’z*diyy’ 


The 1000 cards of weights and heights were used in this way with the digits 
grouped. There resuited (dimensions omitted) 


N = 1,000 

=x = 139,299 Ly = 678,896 
rx’ = 139,520 Sy’ = 679,420 
Sex’ = 19,722,326 yx’ = 94,929,322 


rary’ = 94,848,036 Zyy’ = 461,885,052 


which gives E'[r] = .4957. The ungrouped 4 place value is .4952. 
For use without Hollerith machines, this method indicates the recording of 
the values of yz and x, as well as f., for each entry in the correlation chart. 
The generalization of method 4 leads to 


Lay = U(x’ + e)(y’ + 9) = —2Ta’y’ + Ta'y + Try’ + Len 
so that we have 
(25) E“ [zy] = —Zz2’'y’ + Izy’ + E2’y. 


It is to be noted that the quantity 22’y’ is the unbiased estimate of Uzry 
resulting from the usual frequency distribution. This formula can be used 
with formula (11) of section 4 to obtain an estimate of the correlation coef- 
ficient. 

The correlation chart application of method 4 demands the triple entry 
fey, Ty, yz for each of the squares of the correlation chart. From these 
values it is possible to compute all the entries needed to use method 4. 
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In general the values of E[=2’y’*F,G,] can be worked out with the repeated 
use of lemma 21. The reader who understands the developments of sections 
3 and 4 should have little difficulty in writing out the formulas resulting here. 

It should be pointed out that, in cases where the first and second order mo- 
ments only are desired, it is frequently advisable to avoid grouping by using 
modern computing machines and, in this way, to eliminate the trouble and the 
errors caused by grouping [6]. 


8. Conclusion. There are additional points which might be considered, but 
they would take considerable space and the presentation is now sufficiently 
complete to enable one to obtain some perspective on the proper use of the 
new methods. 

If precision is not needed, the use of the former grouping methods is advised. 
But if additional precision is needed, and if the results of a single grouping 
only are available, it is advised to use the newer methods. Method 2 is much 
more satisfactory than method 1 and, in many cases, will be sufficient, but, if 
additional precision is demanded, one can use method 3 or one of the com- 
pound methods. 

In general there are two kinds of groupings. One is a recorded grouping, 
and expresses the measures in terms of the units which are desired, while the 
second is a computational grouping which is introduced for the purpose of ease 
of computation. ‘Now the recorded grouping, no matter whether obtained from 
discrete or continuous data, is necessarily discrete. Thus the weights, when 
measured, have to be recorded to the nearest pound, or the nearest tenth of a 
pound, or to the nearest hundredth of a pound, ete. The formulas to be applied 
to the results of computational grouping are the formulas for discrete variates. 
If in addition one wishes to correct continuous data for the recorded grouping, 
he may then apply the usual Sheppard’s corrections for continuous data. How- 
ever, it is advised to make the recording grouping sufficiently detailed so that 
the errors are slight. Thus one might record the values of heights to the nearest 
tenth of a pound, but use ten-pound intervals in making calculations. In this 
case the values when corrected for the computational grouping (to the nearest 
tenth of a pound) would presumably be sufficiently precise so that the additional 
grouping for recording would not be necessary. (In many cases the two group- 
ing corrections are combined in a single grouping correction for continuous 
data.) : 

It appears that it is not sufficiently satisfactory to continue to record the 
results of grouping in the usual form of a class mark (or class limits) and a fre- 
quency if the resuits are to be used by others. The table should include an 
additional column of z, and preferably a column of x3, , where the x’ are the 
computational grouped values and the x are the measured values recorded to 
a considerable degree of accuracy. The arrangement takes little more space 
than the present frequency distribution, and it can be obtained from the recorded 
values with a reasonable amount of additional work. In the case of correlation 
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it is suggested that the present grouping of frequencies in the correlation chart 
be augmented with the values of xz, and y, for each square. In this way it 
is possible for those who may use the distributions later to obtain much better 
estimates than would be possible from the frequency distributions as now 
recorded. This point certainly should be considered by all those who prepare 
tables for general use, and yet are forced by practical considerations to use some 
sort of grouping in reporting the results. 
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ON THE CORRECT USE OF BAYES’ FORMULA 


By R. v. MIsEs 
Harvard University 


The problem that we try to solve by using Bayes’ formula consists in making 
an inference from an observed statistical value upon the unknown value of a 
parameter, and in examining the chance of this inference being correct. One 
may call this the principle problem of practical statistics or the estimation 
problem, or, as the author put it in German (Rueckschluss-Wahrscheinlichkeit) 
problem of inference probability; at any rate we encounter this kind of problem 
in various forms in almost every branch of statistical investigation. It will be 
convenient to base the following discussion on a concrete question in quite 
specified form which will allow us to see clearer the points that are to be stressed 
in this paper. 


1. The problem. In examining the quality of water supplies with respect 
to the number of bacterias of a certain kind they contain, a definite procedure is 
usually adopted. One takes n = 5 samples out of the water, each sample of 
exactly 10 cem. Then by a certain biological test one finds out whether or not 
each sample contains at least one bacteria of the kind under consideration. The 
number x (zero to five) of positive tests is the observed value from which an 
inference is drawn upon the probability @ for a sample containing at least one 
bacteria. It is assumed that this 6 is connected with the average number A 
of bacterias per 10 cem by 


(1) @6=1-—e*; 6@=6=063 fori=1 


according to Poisson’s law. <A particular question which we want to answer is 
this: What is the chance of being right, if we- conclude from the observed fact 
xz = O, (in other cases from x = 1) that @ lies between 0 and 6, = 0.63 (or A 
between 0 and 1)? 

For a given @ the probability of getting x positive tests out of n tests is 
according to Bernoulli’s formula 


(2) plz |e) = (")ea — a, 


The chance of having a 6-value between 0 and 6, when z positive tests are ob- 
served is according to Bayes’ formula 


61 
I p(x | 6) dP(@) 


(3) P, (0;) = 


1 


I p(z | 6) dP(@) 
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where P(6) is a distribution function, monotonically increasing from 0 to 1 
and usually known as the a priori probability. 


2. The apriori. The function P(@) is generally considered as a troublemaker. 
As one uses to call P the a priori probability most people think that it has some- 
thing to do with those absurd conceptions of non-empirical, a priori known 
probabilities that cannot be tested by any experiments etc. This cannot be 
strongly, enough refuted. In our particular case the meaning of P(6) is the 
following. Each probability statement refers, as we know, to a certain infinite 
sequence of experiments or trials which form a kollektiv. If we ask for the 
chance P,(6,) of having a 6-value between zero and 6, when a certain x has been 
observed, we have in mind a sequence of trials each consisting of two steps, 
first, picking out one particular water supply, and then testing the number z 
of samples that contain bacillas. Among the first N trials of this kind we shall 
have N, cases where the 6-value for the water supply picked out lies between 0 
and 6,, then we shall have N, cases where the number of positive tests is z, 
and finally in a number N;, of cases both conditions will be fulfilled. The 
chance P,(6@,) we ask for is then by definition 

Niz 


(4) P,(61) _ lim N, ’ 


while the so-called a priori probability is 


(5) P(@,) = tim ™. 


N 
Later on we shall also use the probability 


(6) 


All these magnitudes are to the same extent empirical or non-empirical. They 
are “empirical,” since we get approximate values for them out of a long sequence 
of experiments, and they may be considered as something super-empirical since 
the concepts of an infinite sequence and of a limit are used in the definition—as 
each theory must involve a certain amount of “idealization.” 

In order to avoid the above mentioned equivocation the author had sug- 
gested a long time ago’ to call the probabilities corresponding to P(@) and P.(6) 
respectively the initial and the final probability. Another expression which 
could be used in connection with the distribution function P(6) is overall distri- 
bution, since it means the distribution of @-values within the total mass of 
samples, not regarding what the values of z are in each case. 


3. No randomness required. Now, the first remark we have to make is the 
following: In the Bayes’ formula (3) the existence of a function P(6) is presup- 


1 Cf. reference [2], p. 152. 
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posed, i.e. we assume that in the sequence of successive trials the frequency of 
those cases in which @ falls into a certain region has a definite limit. But 
nothing is assumed about this limit being independent of a place selection. 
The sequence of trials must fulfill the first condition of a kollektiv, with respect 
to @ but not the second; in other words the randomness in the succession of 6-values 
is not required. Thus we may say that @ is not supposed to be a chance variable 
in the usual sense of this term. Sometimes people are shocked by the idea that 
in Bayes’ theory the individual cases are supposed to be picked out at random, 
and it is often considered as a superiority of the method of confidence intervals 
that here such assumption is avoided. 

It is true that in the latter method even the existence of the frequency limit is 
not required,’ but this does not seem to make any essential difference. The 
fact is that, if we want to make an inference upon the value of @ i.e. an assertion 
about the chance of @ falling into a certain interval, we have to assume that in 
the long run different @-values may occur with certain frequencies. 

It may be useful to have different expressions for the two cases where a fre- 
quency limit is or is not supposed to be independent of an arbitrary place 
selection. As we use the word probability in the first case it seems suitable to 
apply the word chance in the second. Thus, if P(@) is the initial or the over all 
chance of 6 we would say that P.(@,) is the final chance of 6 being smaller than 
or equal to 6, for a certain observed z-value. When P(@) is supposed to be a 
probability, i.e. to fulfill the condition of randomness, then P,(6;) will have this 
property too and has to be called probability. 


4. Inequalities for the final chance P,(6). A much better founded objection 
against the practical application of Bayes’ formula consists in saying that in 
most cases we have no sufficient information about the function P(6). This 
undeniable fact leads often to an incorrect simplification of the formula by re- 
placing in it dP(6) by dé which means an a priori probability of constant density. 
It is obvious that this is no solution: if you do not know what P(@) is, to assume 
it equal to 6. On the other hand, if we accept Bayes’ formula as correct (and 
there is no reason for not doing so) we learn that the value P.(6) we ask for 
depends essentially on P(6), and is undetermined as far as P(@) is undetermined. 
The only consequence in this situation is, first to use all information we can get 
about P(6), and then to make the answer as vague or undetermined as the in- 
completeness of this information requires. 

One way to do this consists in setting up inequalities for P,(@) based on 
certain inequalities for P(@). A formula which turns out to be useful, at least 
in a well-known asymptotic problem is the following: 

Let us consider the general case where 6 stands for several variable parameters, 
and let A be the set of all possible values of 6. We are interested in the final 
probability P; of a subset C of A given by 


2Cf. reference [4], p. 201. 
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I |, P(e |6) dP) 
(6) fi ee i cceeemeniice 


[ (, ir) aPC) 


where x is supposed to be known. 
Let P’ be the value of P, under the assumption of a constant initial density 


and denote by Ps , P’; the analogous values for a subset B which includes C so 
as to have 


(7) C<B<A. 


The quantities P, and P, depend only on the function p(x | 6) and the sets B 
and C while Ps and P¢- change with P(6). 
If we assume that the initial density p(@) has the limits 


(8) m < p(0) <M within B 


m’ <p(0) <M’ within A — B, (A minus B) 


it can easily be shown that 


(9 ~Pe+Ma-Ppattst r+ a- Py. 
Cc 

We may consider the following application of these inequalities. 

If we are concerned with a case where a great number n of trials is involved, 
the function p(x | 6)—which determines the P’ values—shows an increasing con- 
centration at a certain point of the set A. In other words, for large n we have 
a subset B more and more reducing to one single point for which P, is as near 
to 1 as we want. If we then assume that the density p(@) is continuous and 
bounded, the difference between m and M tends to zero, and if m is supposed to 
have a positive lower bound, both the first and the last expression in (9) tend 
to unit or P; approaches Pe. This is a generalized form of the statement 
which the author proved for the first time in 1919,* that in the original Bayes’ 
problem where we are concerned with n repetitive observations of an alterna- 
tive, the final probability becomes more and more independent of the initial proba- 
bility P(6) as the number n of observations involved increases. 


5. Using previous experience. The inequalities (9) may be of use in many 
cases. But to be sure, in general, they are not the basis upon which practical 
estimation judgments rest. Everybody acquainted with the conditions of test- 
ing water supplies takes it for granted that the outcome x = 0 (no positive test) 
supplies a sufficient reason for the statement 6 < 6, = 0.63 (less than one 


3Cf. reference [1], p. 84. 
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bacteria per 10 cc). But, if nothing were known about the initial distribution 
P(@), we could assume P(@) in the form 


P(6) = 6", p(@) = me” for0 <6<1, 


with a large value of m. With n = 5,2 = 0 equations (2) and (3) give Po(@,;) = 
0.50 for m = 10, and Po(@:) = 0.88 if mis 5. These values are much too low 
to justify any recommendation of a water supply for which x was found to be 
zero. Thus we have to ask: What is the real source of the confidence we put in 
the inference from x = QO upon @ S @,? 

There is no doubt, that this confidence is based on previous experience. We 
know that the water supplies subjected to the routine test in the past formed a 
class of rather clean than dirty water and we rely that a new sample will belong 
to the same class. The author was given the following information about the 
results under the jurisdiction of Massachusetts during the last decade. Out 
of a total of N = 3420 examinations there were found 


3086 cases with x = 0 (no positive test) 
279 cases with x = 1 (one positive test) 
32 cases with x = 2 
15 cases with x = 3 
5 cases with x = 4 


3 cases with x = 5 


The overwhelming majority of cases with x = 0 is evident. The question is 
only how we can use these statistics of past experiments for obtaining a nu- 
merical inference upon the value of P,(6). 

If the initial distribution P(@) were known, we could find the probability Q. 
of getting x positive tests out of n: 


(10) Q.= [ r@\0 are) = (") [ oa - 9" are. 


Using the numbers N,, N., Niz introduced in section 2 the probability Q(x) 
is defined by equation (6). 

If the number N of past examinations is considered as sufficiently large, we 
can take the ratios 3086/3420, 279/3420 etc. as approximate values for Qo, 
Q: etc. Now, according to the well-known identities 


(11) : > x (") 67(1 — 0)" * = 0. 


n z=0 


(12) a, ee (")era ——"* =f, 


n(n = 1) z=0 
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and using (10) we can derive from the values Qo , Qi, --- , Q, the first and second 
moments of the distribution function P(é@): 


1 1 n 
(13) 


n 


MM, = [ 6° dP(6) = 2 > x(x — 1)Q-. 


n(n 
If we introduce here the above mentioned empirical ratios for Q, we find the 
approximate values for the first and second moments of P(é): 


(13’) M, = 0.02474 M, = 0.00401. 


6. Determination of a distribution function by its first moments. In an 
earlier paper the author showed [3] how the exact upper and lower bounds for a 
distribution function P(6) can be found, if the expected values of two functions 
j(@) and g(@) are known. The only condition was that the curve represented 
in a Cartesian coordinate system by x = f(6), y = g(6) is convex. Let us take 


f(®) = g(6) = 0 for@ <0 
(14) fo =0 g)=6 forO0s6<1 
f(0) = g(@) = 1 for @> 1. 


In this case the condition is fulfilled and the expected values of f(@) and g(6) 
are the moments M,, M2, respectively. The results obtained in the paper 
quoted above take the following form: 

First, we have to derive from the given values 1/, and Jf, two points 6’ and 
6’ of the internal 0 < 6 < 1 








em M, - MM, 0 M 
(15) ies i-W,’ a M,° 
Then the limits for P(@) are: 
M.~ Mi 
te esac ee ae <6@<@' 
05 PO Spore forOZ6<6 
(16) i- M,-— =oS *< Po) S1- oe = ; foré’ <6 < 6” 
_—Oh- 8"  < p@ <1 for 6” <6 <1. 


M. — 2M,0+ 6 ~ 


In our case we find 6’ = 0.0213, 6” = 0.1619 and the point 6, = 0.6321 falls 
into the third interval 6,1. ThelinesO A BC andODE F G in Fig. 1 show 
(slightly distorted) the lower and upper bounds for P(6). 





o.033 
0.002 


} » 2 3 

Fre. 2 

Fig. 1. The limits of the overall distribution function 
Fig. 2. The 99% region in the methods of confidence intervals 


7. Application to Bayes’ formula. The inequalities (16) enable us to find in a 
simple way a lower bound for the end probability P.(6@,) defined by (2) and (3) 
in the case x = 0. Let us denote by A the numerator in (3) and by B the 
supplementary integral 


(17) B= I p(x| 6) dP(6), 


so as to have A + B for the denominator in (3). If the subscripts min and max 
denote a lower and upper bound respectively we can write 


A Amin 
- ee ere ae ee 


Now, taking x = 0 we find by product integration 
61 
(19) A= P(@)(1— 6)" +2 I P(@)(1 — 0)" as. 


Therefore, Amin is found when we introduce in this expression the lower lim 
for P(6) as given in (16). If we do this and use the values for M, and M, 
according to (13’), numerical computation leads to Amin = 0.712. 

In the same way we obtain B in the form 


(20) = —P(,)(1— 6)" +n [ P(0)(1 — 6)" a8. 


The upper bound Binax is reached, if we introduce in the integral P(6) = 1 and 
in the first term the minimum value for P(6,) following from (16). The second 
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term becomes thus equal to (1 — 6)” and the numerical result is Bmax = 
0.0000607. Therefore the inequality (18) supplies 


0.712 
71206 


The final outcome secured in this way can be formulated as follows: If we 
assume that in continuing the experiments the distribution of test results will be 
about the same as it has been in the past 3420 cases, we have a chance of more than 
99:9% of being right, when we state in each case of no positive test that the density 
of bacterias is less than 1 per 10 ccm. 

The high value of 99.9% for P(@,) is of course strictly bound to the assump- 
tion that the entire mass of water supplies to be tested is homogeneous and 
sufficiently characterized by the distribution of test results found in the past. 
If e.g. we had to assume that the six possible values for x (0 to 5) in the long 
run appear with equal frequencies so as to have Qo = Q; = --- Qs = 3, the 
same method would give M, = 3, M2 = }, then 6’ = 3, 6” = #, and the final 
result would be Po(@:) 2 0.73. The assumption of a constant initial density 


P(6) = @ would give Po(@:) = P,(6:) = 0.9975, a little less than the value 
found above in (18’). 


(18) Po(@1) = 0 = 0.99915. 


8. The case x = 1. The results are less favorable in the case of one positive 
test, x = 1. Here we have 


(21) p(1 | 6) = neil — 6)” = 50(1 — 8)", 


and the derivative of p is first positive, then negative. We can conclude from 
Fig. 1 that the minimum value for A and the maximum for B will be reached 
when the distribution function P(@) is represented by the line O DI H JG 
where J H is horizontal and H the point on B C with abscissa 6,. The abscissa 
6) of I is determined by the equation 


(22) — Mi— Mi _ (i — 6) __ 
M2 — 2M, + @ M, — 2M,6, + 62’ 


which supplies 6 = 0.0190. We then have 


8 
(23) hein * I p(1| 6) aP(6), 


with the value p(1 | @) from (21) and with 


— M—-M 
PO) = Mz — 2M,0 + @ 


according to (16). On the other hand Bmax is found, as in the former case, to be 


(24) | = p(l | 6,){1 ae P(6)], 
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where we have to take for P(6,) its minimum value according to (16). The 
numerical computation yields Ajnin = 0.0062 and Bnax = 0.00052 so as to give 


> 82 _ Q99. 


The result is that—under the assumption above mentioned—we have more 
than 92% chance of being right, if we predict each time one out of five tests has 
been positive that the density of bacilli is less than 1 per 10 cem.—The chance 
computed under the assumption of a uniform initial distribution P(é) = 6 
would be 0.97. 


9. The method of confidence intervals. One may ask what kind of answer 
to our questions can be deduced from the principle of confidence intervals. 
This method has undeniably to its credit that no use is made here of the initial 
distribution P(6) and that, therefore, all its statements are completely inde- 
pendent of what is assumed about P(@). 

In order to apply this method‘ we have to select for a given degree of confi- 
dence, say a = 0.99, a region of acceptance, i.e. an area in the two dimensional 
x, 6 plane limited by two lines x:(@) and 2x-(@) so as to have for each 6 


(25) Prob {x,(6) S x S 22(6)} = a. 


The region is, of course, not uniquely determined by (25). In our case, how- 
ever, one will generally agree that the best way to determine the region consists 
in assuming for x:(@) and 22(6) two step lines with steps at the integer values 
x = 0,1, 2, --- as indicated in Fig. 2. Then the formula (2) for p(x | 6) com- 
bined with (25) supplies the abscissae of the steps, if zis given. If we transform 
the limits for @ into limits for \ using equation (1), the final outcome reads as 
follows: 

Whatever the initial distribution P(6) may be, we have a chance of 99% of being 
right, uf we predict: 


each time x = O is observed that d lies between 0 and 0.92, 
each time x = 1 1s observed that d lies between 0.002 and 1.51, 
each time x = 2 is observed that d lies between 0.036 and 2.24, 
each time x = 3 is observed that d lies between 0.112 and 3.41, 
each time x = 4 is observed that X lies between 0.25 and 8.48, 
each time x = 5 is observed that d lies, between 0.51 and ~. 


It is true that in this way we obtain a result independent of any assumption 
on P(@). But it is essential that the chance of a = 99% holds only for the siz 
joint statements as a whole. This means it may happen that for instance the 
first assertion (that A is smaller than 0.92 in the case x = 0) is correct but very 
seldom or even never, while other assertions (e.g. those for x = 4 and 5) have 


‘Cf. reference [5] and reference [4], p. 203. 
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a much greater chance than 99% of being correct. Whether this happens or 
not depends on the initial distribution P(6). As long as we know nothing about 
P(6) we are not in the position to conclude, by using the method of confidence 
intervals, that the particular statement ‘A < 0.92 if z = 0” has a chance of 
99% or even any chance at all of being correct. On the other hand, when xz = 0 
has been observed we are in no way interested in consequences that may be 
drawn in the case x = 4 or x = 5 or in a set of statements that includes the 
cases x = 4 and x = 5. The only practical question that is relevant to the 
purpose for which the tests are made is this: What can we conclude from the fact 
that in a certain instance x = 0 has been observed (or in another instance x = 1)? 
It seems that the method of confidence intervals, discarding any consideration 
of the initial distribution, can supply no contribution towards the answering 
this particular question. 


REFERENCES 

[1] R. v. Mises, ‘‘Fundamentalsaetze der Wahrscheinlichkeitsrechnung,” Math. Zeit., 
Vol. 4 (1919), pp. 1-97. 

[2] R. v. Mises, Wahrscheinlichkeitsrechnung und ihre Anwendung in der Statistik und 
Theoretischen Physik, Leipzig und Wien, 1931. 

[3] R. v. Misgs, ‘‘The limits of a distribution function if two expected values are given,” 
Annals of Math. Stat., Vol. 10 (1939), pp. 99-104. 

[4] R. v. Misgs, ‘‘On the foundation of probability and statistics,’’ Annals of Math. Stat. 
Vol. 12 (1941), pp. 191-205. 

[5] J. Nevman, Roy. Stat. Soc. Jour., Vol. 97 (1934), pp. 590-592. 




































AN ITERATIVE METHOD OF ADJUSTING SAMPLE FREQUENCY TABLES 
WHEN EXPECTED MARGINAL TOTALS ARE KNOWN 


By Freperick F. STEPHAN 


Cornell University and U. S. Census Bureau 


1. Introduction. In a previous paper by W. Edwards Deming and the 
author [1] the method of least squares was applied to the adjustment of sample 
frequency tables for which the expected values of the marginal totals are known. 
From observations on a sample the frequencies n,; for the cell in the 7th row 
and jth column of a two dimensional table and the r row and s column totals, 
n;.andn.;,are obtained. These frequencies are subject to the errors of random 
sampling and it is desired to adjust them so that the row and column totals 
will agree with their expected values, m;, and m.;, which are known. The 
adjustment involves the solution of the r + s — 1 normal equations 


mids. + > nyd.5 = Mis — %., $= 1,2, ---,f 
7 
(1) 


De msds. + 1.52.; ms Mas 3=1,2,---,s—1 

where the \ are Lagrange multipliers from which are calculated the adjusted 
frequencies 

(2) ms; = nj(1 + Ax. + X.)). 

Similar equations arise in the three dimensional case. 

A method of iterative proportions was presented for effecting the adjustments 
more conveniently than by solving the normal and condition equations, and it 
was stated that ‘the final results coincide with the least squares solution.” 
This statement is incorrect, for although the adjusted values satisfy the condition 
equations, they do not satisfy the normal equations and hence they provide 
only an approximation to the solution. The method of iterative proportions has 
several interesting characteristics that will be discussed in a later section. 
This paper now presents a method that converges to the values given by the 
least squares adjustment and is self correcting. It can be used with any set of 
data and weights for which a least squares solution exists. The two-dimensional 
case will be considered first. 





2. The two-dimensional case; expected row and column totals known. 
Assume that a sample of n items is drawn at random and cross-classified in a 
table of rrows and scolumns. As in the previous paper, let n;; be the frequency 
in the 7th row and jth column of the two-way frequency distribution. Indicate 
summation by substituting a dot for the letter over which the summation is to 
be performed. Then n;. and n.; are the marginal totals for the ith row and 
jth column respectively. Let m;, and m.; be the expected values of these 
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marginal totals calculated from other information or from theoretical considera- 
tions, and c;; a set of constants known or estimated to be proportional to the 
reciprocals of the weights of the n;;, i.e. proportional to their error variances. 
Since the weights are positive, the c;; are non-negative and finite. It is assumed 
that the set of weights is such that for the given data an adjustment exists. 
The least squares adjusted frequencies m;; can be computed from the given 
numbers ¢;;, n:;, m;., and m.; by a series of approximate adjustments in a 
manner now to be explained. Let m{?) be the pth approximation to m;;. In 
conformity with this notation mS) = .n;;. Let 
(3). dP =m ;— mi, di? =m.—m?, dd? =m; — m?, 
be corrections that must be added to the m‘”) to produce the least dquares 
adjusted frequencies. As d — 0, m®” — m. Let d{”) and d°%? be constants 
determined arbitrarily between the limits set by equations (5) to (7). Any one 
\ may be fixed arbitrarily and kept constant through successive approximations. 
Note that AS” = A = 0 and that, if at every step we set A“? = 0, the \” 
are approximations to the Lagrange multipliers in the normal equations. After 
p steps in the iterative process the approximate adjusted frequencies will be 


(4) ms? = nis + c1(M” + AY). 
The following conditions, derived from (19), (23), and (24), are sufficient to 


make the successive approximations converge to the least squares adjusted 
frequencies: 


val = 
ri” = nf? , + o§?) ds? Tex. ? 
p p—1 p—l 
ry ) = ng ) + eo” ad‘ Ie; 


(6) o<a”, ose, arr+o? <2, 


(5) 


and, for at least one pair 2, 
(7) AP (ary + Pay Y>0; a + a” <2. 


The 6’s are introduced because in actual computations the successive approxi- 
mations \‘” can only bé calculated to a limited number of digits and because 
the adjustment may progress more rapidly if the computer is permitted to use 
his judgment in determining the approximations as he observes the ‘course of 
previous approximations. 

The process of adjustment is continued until the d{”’ and d‘?’ becorfe suffi- 
ciently small to provide the desired degree of agreement between the adjusted 
and expected row and column totals. 


3. Example. The following example shows the steps in the adjustment for 
a table of 3 rows and 4 columns with 6” = 0°”. = 1: 
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i | is iia | Pi | ds | ofa | a” | a? | a yes | (2) n® a? ; a 
| ; 4) 8) | © | | (10) 
un) mi —|—| wm—| me—| — | — | mae — | m 
12| 7426) — | — | 455 — | 7505.6 — |; — — 7496.3} — | 7497 
13| 4709) — | — | 358} — | 47126) — | — | — 4709.6 — | 4711 
14} 2145) — | — | 176, — | 2055.83) — | — | — 2051.3} — | 2049 
| | | | | | | | 
21; 517 — | — | 52 — | 52899 — | — | — 529.5, — | 529 
22| 928) — | — | 95) — | 973.47 — | — | — 978.3} — | 979 
= oe i - |) we] ae Kt Kl 643.1} — | 644 
4) 703) — | — | Wi — | 688.7; — —- | - 691.9} — | 692 
Doe 4 | 4 | } | 
ai mi-i- | wa | Oe | CU _ 201.1) — | 201 
32| 373) — | — | 381 — | 369.1; — | — | — 372.3} — | 373 
$3); 337} — | — | 3] — | 32.7, -— | — — | 831.7) — | 332 
4) 4235 — | — | 3 — 394.5) — — — | 397.5) — | 397 
| | | | | 
| } { | | | | 
1) 1507| 1501) 6) 146/—.041; 1506.7; —5.7| —.0390| —.0800| 1503.5, —2.5, 1501 
.2| 8727] 8849| +122) 588/+-.208) 8848.1) +0.9| +.0015| +.2095) 8846.9) +2.1) 8849 
.3| 5668) 5687) +19| 445/+.043) 5680.8) +6.2) +.0139] + .0569| 5684.4) +2.6) 5687 
-4| 3273 3138) —135, 285|—.474) 3139.0, —1.0| —.0035| —.4775| 3140.7, —2.7) 3138 
| | | | | | | | | | 
1. | 15063] 15028) —35| 1064 — .033| 15051.5|—23.5) —.0221) —.0551) 15030.1) —2.1)15028 
2.| 2770) 2844 +74 273 +.27) 2830.5/+13.5) +.0495) +.3195 2842.8) +1.2) 2844 
3. | 1342) 1303) —39) 127) —.31) 1292.6/+10.4| +.0819| —.2281) 1302.6, +0.4) 1303 
Ed | | | 
19175 19175, 0) «1464, — | 19174.6, +0.4, — | — | 19175.5) —0.519175 



















Columns (1), (2) and (4) are given. Columns (3) and (6) to (11) are calcu- 
lated in succession using equations (3), (4), and (5). It is not necessary in 
practice to record the 6’s or even determine their values since the \‘” may be 
determined directly at convenient values approximately equal to their corre- 
sponding d{?~” + d§?~”/c;. and AS?~” + d‘?~-”/c.;. The final adjusted fre- 
quencies given in column (12) are derived by another repetition of the adjust- 
ment process but the amounts involved are so small that they can be calculated 
mentally and the n;; rounded at the same time. 














4. Computing procedure. The computing procedure may be set up in any 
of a number of ways to meet the preferences of the computer and the charac- 
teristics of the problem. Ordinarily it is desirable to make every number 
positive and the procedure as nearly routine as possible. 

For two-dimensional adjustments the following procedure of computing alter- 
nately by columns and by rows is convenient: 

(a) Set up a table of the c;; in r rows and s columns. Enter the c;, in the 
s + 1 column, the c,; in the r + 1 row, and c., = >> c;, = )>c.; in the com- 

. 2 


mon cell. 
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(b) Calculate the quantities A;, = (d{/c;) + a and A.; = (d%/c,)) + 
and enter them in the s + 2 column andr + 2row. The constant a is selected 
at some value that will make all quantities in the computations positive and 
may be any convenient integer greater than 2 max | d{°”/c;. | or 2 max | d?/c,; |. 

(c) Calculate the factors u{” approximately equal to the A;. — 3a and enter 
each on its corresponding row in the s + 3 column. Throughout the computa- 
tions the u»” are merely ),; + 43a. 

(d) Take column j and multiply each c,;; by its corresponding y{” accumu- 
lating the products in the calculating machine. Divide the sum of products 
by c.;, subtract the quotient from A.;, and record the difference y in the jth 
column on the r + 3 row. Repeat for each of the other columns. 

(e) Take row i and multiply each c;; by its corresponding y? accumulating 
the products in the calculating machine. Divide the sum of products by c¢;,. , 
subtract the quotient from A;,, and record the difference »{° on the ith row 
in the s + 4 column bordering the table on the right. Repeat for each of the 
other rows. 

(f) Repeat steps (d) and (e) alternately until a satisfactory degree of stability 
is reached in the u;. and u,;. Then compute each adjusted frequency as follows: 


(8) ms? = cx s(us? + uP —a)+n;, 


taking either u{?? = u{?~” or nS? = u‘?~” as the case may be. 
(g) The computations may be checked at any step by computing 


(9) Lewy? 4 = DP Aies— Do wl? a. = ac. — Dwi? cz., 
7 ‘ ‘ 
or 


(10) dX wi? c, = Z Ai.cG. — d a Cj = a,, — pm wp? Cj. 


i 7 


(h) At any step a constant may be added to all the u{” and subtracted from 
all the ‘?; this may be necessary to keep the y’s all positive. It has no effect 
on the value of a to be used in (8). 


(i) If it is desired to “inflate” the adjusted frequencies (2 mis # 2d Ni 3) 
first multiply each n;;,;. , and n,; by the factor 2d mM; Id, Nj and then gencced 


as above using the products in place of their mauling Nj, 2. and n.;. 
(j) If before the iterative process has reached an acceptable adjustment it is 
desired to force a satisfaction of the condition equations, compute: 


(11) m2t = ous? + uP — a) + riz + (iP; + dPex.)/c.., 


in which either the d{”’ or the d‘?? are all zero. 


5. Adjustments in three dimensions. If the sample is cross-tabulated in a 
three-way frequency distribution, there are two cases that are not reducible to 
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two-way distributions. These are designated Case III and Case VII in the 
earlier paper [1]. The adjustment equations are, respectively, 


(p) (p) ) (p) 
mie = Nise + ci(As? + AS? + ASP?) 

(p) (p) (p) (p) 
Mire = Nije t cin(AG? +A + A7’), 


subject to conditions on the choice of the \ corresponding to equations (5), (6), 
and (7). For Case III, the conditions are that 


(13) 0< &, o< oF, O0< of, of + of + of <2, 


and for at least one triple ijk, 6? (d{?-”)? + oP (de)? + oP (a?) > 
0 and 0?) + 0%? + o'?? <2. Similar conditions apply to Case VII. 

The computing procedure described in Section 4 can be extended readily to 
the three-dimensional case. For example, in Case VII calculate »$}? approxi- 
mately equal to (d$%) /c;;.) + 4a and u{? approximately equal to (d§& /ci..) + 4a. 
Then multiply each c;; in the column jk by its corresponding (u$} + p!? 
accumulating the products in the calculating machine. Divide the sum of the 
products by c. and subtract the quotient from (d‘/c.;.) + a. Record the 
difference as u‘ and repeat the process for every other.jk column. Take 
us = us} and repeat for each ik column to obtain u$? ; then take un = uQ 
and repeat for each 77 column to obtain us; and so on. The final adjusted 
frequencies are 


( 
(14) mie = nije + cize(us? + ul? + uw? — a). 


(12) 


6. The general case. The iterative method can be extended readily to 
more than three dimensions and to various systems of condition equations. A 
simple general notation may now be introduced. Let the cells be numbered in 
any order from 1 to ¢ and for the 7th cell let n; be the value given by the sample, 
c; a finite positive constant known or estimated to be inversely proportional to 
the weight of n;, m; the least squares adjusted value to be determined, m{” 
the pth approximation to m; , d‘” = m; — m{”, and mf” = n;. Assume that 
the values m, of certain linear combinations of the m, are given, i.e. there is a 
system of consistent linear equations of condition numbered in any order, the 
oth equation being 


*(15) Z biem; = Me ; » bie > 0, 


b;, and m, being known a priori. The corresponding linear combinations of the 
n; and d{” define 


(16) Ne = DLibiem, de” = Di died {”. 
Let 
(17) Co = Dy viet. 
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The pth approximation to m; is 


(18) mi” = ni + c Do bids”, 
where 
(19) = PY + Ody Y/ee, dre” = 0, 


the 6”, and therefore the d{”, being arbitrary for a finite number of steps, 
say p’, but determined thereafter so that 


(20) 2 >> 0? (d?-?)*/ee — Do eld. die 06” dl?” Jee)? > (d’”-”)?/(c, H), 
e q e 


r being a value of c, chosen at the pth step, for which (d{?~”)?/c, is a maximum 
and H a finite number greater than 1 fixed prior to the first step as large as one 
will. That this condition can be satisfied may be shown by putting 0!” = 
1 and 0” = 0 (o # 7). 

A weighted average of several of the possible selections of 6{” satisfying (20) 
will also satisfy (20), positive ‘“‘weights” being assumed. Let k added to the 
superscript represent the kth such selection and let a” > 0 be a constant for 
“weighting” the kth selection in the weighted average which may be chosen 
arbitrarily except that » a?” — 1, Then, if the kth selection of 6{” is repre- 


sented by 6S?”, the weighted averages are 057° = >> a‘? 6”. Substitute 
k 


them in the left-hand side of (20), 
2 lee — Des (2, bel?” de? Jen)? 
@) = 2D LaMar (alryt/e,— De LDbval™ aay”) 
- » a?) (2 > gfe (a>)*/e.) _ > o( 2 al?) Zz: bie afr gir") /e,)?, 


which by the Cauchy-Schwarz inequality 
> (p,k) 2 ef?» = 2 . 
> Dal (2Y arr (ay”)*/e) 


ad pe c( >. ¢*) { Ya” (>> bie afr gir 7¢,)*} 
+ k k ¢ 

ia Z. al?” {2 7 gf? (gi?-?)? Je, a Zz. C; (> Die af?®g’?- /c,)*} 
k oe s o 


2 Deal (dy? )*/ (eH) = (d;”-”)’/(c,H). 


A simpler and more convenient but somewhat more restrictive condition may 
be derived as a special case of (20). Let 0S” = 0 except for a set of one or 


ee 
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more o so selected that b,.-b:., = 0 for every 7 and every pair o’ and o”’ in the 
set. Then (20) becomes 


(22) > {206 — (aS”)?} (dS?-”)?/e, > (d’?-”)*/(c, H). 









Differentiating partially for a maximum with respect to one of the 6S”, we find 
that this special case of the condition will be satisfied if for one o in the set, 
say 7, such that 


(23) (dy? )?/ee > (dy? -”)*/(cr/ A), 
the value of 6°” is chosen in the range, 
(24) 1/(2VH) < 6,” < 2 — 1/(2VH) 





















and for every other o in the set 
(25) 0< a” <2, 


all o$” not in the set being zero. A weighted average of such values of @ will 
satisfy (20) whence (6) and (7) follow. 

In practice values of 6S” satisfying (20) may be selected conveniently by the 
following procedure: 

(a) Select a set of o for at least one of which 6 satisfies (23) and for every 
pair of which b;.-b;,.., = 0. In so far as this restriction permits choose the ¢ 
corresponding to the larger values of (d‘?~”)*/c, . 

(b) Determine values for each 0” in the set approximately equal to 1. 
Until other values are assigned to them assume all other 0$” = 0. 

(c) Choose a o not in the set, say p, for which (d‘?~”)*/c, is relatively large 
and select a value for 6S” such that 


(26) a” = {ar — De Cidip die Beda” /¢e} /d,”™”. 
t op 
(d) Having changed 6$” from 0 to a value approximately satisfying (26), 


continue with other o not in the set letting p in (26) represent each in turn. 
The work may be terminated at any stage leaving some 0{” = 0. 












7. Convergence of the adjustment. 
in the following form 


(27) dD. beds” = d?”, 


The condition equations may be written 


as a system of consistent, but not necessarily independent, linear equations. 
They may also be written as conditions on the m;. The least squares adjust- 
ment minimizes the quadratic form 


(28) Ss? = bo a)*/e 
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subject to the restraints (27). Since the c; are positive, S“ is pesitive definite, 
and therefore a minimum exists and is non-negative. The values of the d{” 
that minimize S while satisfying (27) are m; — n;, the n; being known and 
the m; being the least squares adjusted values that are to be calculated. 

If r is the rank of the matrix || b;. ||, then from (15) and (16) it follows that 
r of the d<” may be expressed as linear functions of the t — r other dS. The 
latter then constitute a set of t — r independent variables. The normal 
equations 


(29) as /ady” = 0, 


are obtained by differentiating S“ with respect to each one of them in turn, 
one equation resulting for each value of h corresponding to a d; in the set of 
independent variables. The normal equations (29) are a system of t — r 
independent linear equations and can be written in the form 


(30) Landy” = TB”, 


where the first summation is over the set of independent variables, and the 
second over the d{°” in the r selected condition equations. The right-hand 
terms are constants. Since a least squares adjustment exists the equations are 
consistent and the rank of the matrix || esa) || is ¢ — r. Any d{° in the set, 
say dS’, is the quotient of two determinants the divisor being the determinant 
|aiay | and the dividend being the determinant obtained by replacing the 
ain) by >> Bea ds”. Consequently each d\” whether in the set or not is a 


linear combination of the d°° and the sum of the absolute values of the coeffi- 


cients of the d°” is finite. Therefore 
(31) max | dj”/~/e; | < G max | d?/Vc | 


where G is (max c,/min c;)' times the sum of the absolute values of the coef- 


ficients of the d® in the linear combination for which such sum is a maximum. 
From (28) 


(32) S® < t max {(d{)*/c;} < G’t max {(d\”)*/ce} 
whence 
(33) (d\”)*/ce, = S/(G’E). 


Consider now the discrepancies 


1 ) —1) 
(34) ds” = m — mi” = d?— — c >) die 8” ds” Jee 
go 


between the m; and the corresponding approximations m{” and write the 
quadratic form 


(35) Ss” = Do (di”)*/e. 
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From (16), (18), and (34) 


(36) di = d +e; >) bie dd”, 
and 
(37) ds” = ds” + DoD biebiycidy”. 


Hence the substitution of (36) in (27) merely changes (0) to (p) in the super- 
scripts, the new equations being consistent by definition and the corresponding r 
of the d{”’ being expressible as linear functions of the other t — r. Further 
(35) is positive definite and hence has a minimum, in fact substituting (36) in 
(28) we find that 

3 so a so a 


so = arte = aw IS +22. Didi? beds” + Docc (2D biede”)*} 


ade ~=—s adk””— ss ad” 


_ @ 
i ad” 


(38) (p) 
(s +2 > a” {?) a as’ as 
a - o o adi” 


Hence a least squares solution for the d{” exists and it leads. by (34) to the same 
values for the m; as does the solution for the d{”. Since the coefficients aia) 
and 8.x) and the number G are functions of the b;, and c; they are invariant 
for the substitution. Consequently (30), (31), (32), and (33) may also be 
written with (p) in place of (0) in the superscripts (33) becoming 


(39) (d;”)*/ce, > S”/(G't). 
From (20), (34), and (35) 
Ss? =D (a\”)*/e; 
= dir )*/es — 2 De di? die a” de?” /Cy 
(40) + Docs (DX bie 96” ds?” /ee)’ 
= SP — 2 oP d)*/00 + Dc (DL bua li'ds-” fon)? 


< s?-? — (d?-”)’*/(c,H), p>p' 
and from (39) 


(41) Ss” < sev is S?-) /M < st ae 1/M}?-”" 
where 
(42) M = G@H/t. 


Therefore, as p > ~, p — p’ > ~, S’” — 0, di” — 0, m{” — m, and conse- 
quently the successive adjusted frequencies obtained by an iterative process in 


— 
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which condition (20) is satisfied converge to the adjusted frequencies that are 
obtained by solving the normal equations. 


8. Rate of convergence. The computer is not as much interested in the 
proof of convergence as he is in how rapidly the successive adjustments reach a 
satisfactory degree of approximation. Equations (39) or (41) are of no help 
to him. The adjustment may be made in one step, with every @ = 1, (a) if 
the condition equations are such that every bib;... = 0, 0’ ¥ o”, ie. if the 
adjustment can be separated into one-dimensional cases when redundant condi- 
tion equations are ignored, or (b), in the two and three-dimensional cases, if 
the c;; or c;;, are proportional to the c;. and c,; or to the ¢;.. , ¢.;., ¢... OF Ci;., 
Ci. , and c. , respectively. Except in these and other special cases the rapidity 
of convergence depends on the d<° as well as on the || bic; || matrix. However, 
it seems that one can make very little use of the d{” to determine the rapidity 
of convergence without actually computing the successive adjustments or making 
some equivalent calculation. 

Certain results can be obtained from the || b;.c; || matrix alone. Returning 


to the two-dimensional case and its notation, consider the matrix || ¢;; || and 
define 


(43) bij = Cj — C.€.;/¢..,¢.. = » 6.4. 
d 


Let the adjustments be made with the restriction that 6!” = 0 and @{” = 1 
when p is even, and 6{”? = 1 and 6°” = O when pis odd. Then.if p > 1 


ds? = —D) (ei;/e.)dP-” = DY Dd (eii/e.)) (cii/er.) dy? 
j , 2 


_ = 2X 2X (8:;/c.;) (54;/ey.) ar (f = 1,2, ---,7) 
The sum of the absolute values is 

(45) dil di?| <br Lai? | < or De di? | 

where 

(46) bi= 2 Dal bi/es| v5 

.; being the greatest of the | 6;;/c;. | in the jth column. Similarly for p > 2 
(47) 2d | dP | <b Dd? | soe De | a | 

where 

(48) ba = Ds Do] bii/es. | v5. 


y;:. being the greatest of the | 6;;/c.;| in the 7th column. 
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Assume again the conditions just preceding (44). Let u;, be the minimum 
ex;/c.;in the 7th row. Likewise let v,; be the minimum ¢; ;/c;. in the jth column. 
Then since 2d”) = xd‘? = 0, 


(49) Z| di? | = 22%aP = — 22a", 


the + and — signs indicating that the last two summations are over positive 
and negative values of d‘?’ respectively. When p is even, of course, all values 


of a? = @. 
From (44) 
(50) dj? = > ed? /e.5 = Do cg | dP |/e.5 — Der es | dP |/e, 
? 2 
= om Ci; | in \/c.; —2 _ Cij | aj?” \/c.; 
7 
= 220 ai dy \/es — / ey | dP” |/c.; 
P 2 
and by (49) 
(51) [di |< Docu | dP |/e; — uw. Da?” |, 
7 7 


(52) Lid? | s Dla?” | — dw). 


Similarly 
(53) Dla | <2) air? | (d — doo) 
7 + 7 
Let bs = 1 — De ui. and bh = 1 — > .5, then 
% 2 
(54) Del as? | < bsbe Do | di? | < (baba)? DU | as? |. 


Now 6b; or be may be greater or less than b, or b, but, unlike b; and b2 , they 
can not exceed unity. Let b° be the lesser of b; and b3b,. Then under the 
conditions stated with equation (44) 


(55) Z| dP*? | < S[ dl? | < zl di? | <b” * Z| di | < bP * =| d} |. 
It follows from (40) that 
S” = SP 4 Do di? )?/e.. + De (a”)*/e.3 
‘ 7 


= (Daa. + D@PVes) 


< > (Qu | dS |)?/min c;, + ( | a3? |)?/min c.;} 


< (| ds? | + | af?” |)? {(1/min c;,.) + (1/min ¢,;)}/(1 — v4). 
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The reduction in S‘” in g steps of the iterative process is 


ptq—l 
D= s” = Sirro si 2 [> (d\”)*/e;. + ie (d‘)?/e. 5] 
=p 7 7 


(57) pt+oq—1 
> 2 (CL di? |)*/@ max o.) + QU | di? |)*/(s max e,)]. 


from which, by (55), if g > 1 is odd, 


ui 
sy) D2 WF har | + ar p(t + td. 


r max (Cj. § Max C,; 


The relative decrease in S” is, therefore, by (56), 


D D 1/min c;, + 1/min c,; - 
om Bs Bee> 1+ eae 
S® D+ Ste 
+ b! (b~*? — 1) ( 1 ot ay 
rmaxc; $§ Maxc,; 


If the g steps actually have been taken a better lower limit for the relative 
decrease in S‘” may be obtained by computing D from (57) and using (56) 
for S?*®. Similar equations can be written using by. 

These results can be shown to be valid for an adjustment in which @{?) = 
6‘? = 1 at the first and any of the subsequent steps. They also can be ex- 
tended to the three-dimensional cases but not to three-dimensional adjustments 
with every @ = 1. 


9. Improvement resulting from the adjustment. The least squares adjust- 
ment eliminates a portion of the errors of sampling, i.e. a portion of x°, from 
the set of frequencies estimated from the sample. In fact any adjustment that 
satisfies the condition equations does this. 

Let ¢; be the error in the ith value given by the sample and 6$” the error in 
the pth approximation to the least squares adjusted value. Then 


(60) | 5” =e tc: De bieds”, 
and 
) D (6? ?/eiv = DV /ev + 2D” HH” — Dei (Dvds?) 


The complete adjustment makes 8/?) vanish and therefore, since the last term is 
non-negative, 28;/c; < Ze;/c; except in the trivial case in which all dS? = 0. 
From (37) 


(62) De (65")"/er = De e/eit Do re (ds” — do”). 
The last term may be computed readily at any stage in the iteration. {f the 


sampling is at random, k Y</c; is distributed approximately as x’ with t — 1 
degrees of freedom, where / is the ratio of the c; to the corresponding error 
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variances of the n;. Therefore it would seem appropriate to compute 
k =d, ad, the reduction in x’, as a measure of the improvement achieved in 
the final adjustment. 


10. The method of iterative proportions. The iterative proportions method 
described in the earlier paper [1] implicitly defines, in the two dimensional case, 


(63) Mig = Mi. M.jNij 


the u;. and y.; being given by the r + s condition equations 


(64) ~~ = pm Mi. Mj Nij , a Z Mi. M7 Nij , 
2 t 


any r + s — 1 of which constitute a consistent system of independent equations 
in r + s unknowns. One multiplier, say 4. , may be fixed arbitrarily. Then 
for a2 X s table it is necessary to solve an equation of the sth degree. If s = 2, 
there is only one acceptable solution, given by the positive root; if s = 3, there 
is only one solution of the cubic for which all the adjusted frequencies are non- 
negative. For 3 X 3 and larger tables the adjustment appears to involve the 
solution of equations of the tenth or higher degree and there is then no choice 
but to use methods of approximation. 

The adjusted frequencies given by the method of iterative proportions are not 
identical to those given by the method of least squares. When the adjust- 
ments are small relative to the frequencies adjusted, however, the results given 
by this method approximate those of least squares. For the two-dimensional 
case the successive adjustments converge to a set of frequencies that satisfy the 
condition equations. The author has not found a proof of convergence or 
divergence for more than two dimensions. 


I wish to express my appreciation of many stimulating conversations with 
Dr. W. Edwards Deming on this and related problems, and of the helpful 
critical reading of certain portions of the manuscript by Dr. Joseph F. Daly. 
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ESTIMATION OF VOLUME IN TIMBER STANDS BY STRIP SAMPLING 


By A. A. Hase, 
California Forest and Range Experiment Station’ 


1. Introduction. The present paper is the second of a proposed. series, in 
which it is intended to present a systematic study of the properties of several 
methods of sampling timber stands and statistical treatments of the samples. 

The effects of size, shape, and arrangement of sampling units on the accuracy 
of sample estimates of timber stand volume were reported in the earlier paper [1] 
for 5,760 acres of the Blacks Mountain Experimental Forest. With complete 
inventory data, the nature of stand variation was shown to be such that 2.5-acre 
plots, the smallest size tested, were more efficient sampling units than larger 
plots, i.e., for a given intensity of sampling the sampling error was smaller. 
Long, narrow plots were more efficient than square plots of the same size. 
Line-plot sampling units consisting of two or more equally spaced plots along 
lines of fixed length were as efficient as single-plot sampling units and more 
efficient than strips consisting of plots contiguous end to end. Improvement in 
the accuracy of estimates was obtained by subdividing the area into rectangular 
blocks of equal size, and sampling each block to the same intensity. By sys- 
tematic sampling, whereby the center lines of parallel line-plot or strip sampling 
units were spaced equidistant, the sample estimates of stand volume were im- 
proved over estimates from comparable random samples. Treatment of the 
volumes on individual plots of systematic samples as random sampling observa- 
tions, however, as is sometimes done in practice, was shown to give seriously 
biased estimates of sampling error. 

In the present paper we shall be concerned with sample estimates from strip 
samples taken within blocks of irregular shape, and consequently with sampling 
units which vary in length within samples. The methods will be equally 
applicable to line plot samples. 

Following the general ideas expressed by Neyman [2] it is felt that: (1) If the 
formulae of the theory of probability have to be applied at all to the treatment 
of samples, the theoretical model of sampling must involve some element of 
randomness. (2) This element of randomness may conveniently be introduced 
by a random selection of the sample, but may also be assumed present in the 
distribution of deviations of timber stand volumes in the area sampled from a 
postulated pattern. (3) Many attempts to treat systematic arrangements 
statistically are faulty because the treatment consists in applying to systematic 
arrangements formulae that are deduced under the assumption of randomness. 
If the arrangement of sampling is a systematic one, and random errors are 


1 Maintained by the U. S. Department of Agriculture at Berkeley, in cooperation with 
the University of California. 
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ascribed to Nature, then the treatment of the data should be based on formulae 
deduced under explicit assumption of the svstematic arrangement of sampling 
and of some random element in the material. An example of this kind of treat- 
ment is provided by Neyman’s method of parabolic curves [2] devised for the 
treatment of systematically arranged agricultural experiments. (4) Lastly, a 
mathematical treatment of any practical problem is useful only if the predictions 
of the theory are in satisfactory agreement with the empirical facts. Whether 
the method of sampling is random or systematic, the mathematical theory of 
sampling always involves certain elements that are postulated, either in respect 
to the method of sampling itself or in respect to the material sampled. To have 
a reasonable certainty that a particular mathematical treatment is useful in 
practice it is necessary to make empirical studies to find out whether the devia- 
tions from postulates of the theory that may occur in actual situations do or do 
not seriously affect the validity of the predictions. 


2. Notation and definitions. Before proceeding to the main subject of this 
paper it may be useful to explain the meaning of certain statistical terms and 
symbols. Following Neyman, a sharp distinction is made between three differ- 
ent conceptions that are frequently confused by the practical statistician. 

DEFINITION 1: If wu; , we, +--+ , Uw are any fixed numbers, whether provided 
by some already completed experiment involving randomness, or just arbi- 
trarily selected, Karl Pearson’s term “standard deviation” of these numbers 
and the letter S will be used to denote the expression S = +/3(u; — a)?/N 
in which @ = Lu,/N is the mean of the w’s. 

Now let X denote a random variable, that is a variable the value of which is 
going to be determined by a chance experiment. Thus X may be the timber 
volume on a strip that is going to be selected at random from an area. Denote 
by E(X) the mathematical expectation of variable X capable of possessing values 
Uj,U2,°*:,Un. Then 


E(X) = upr + uspe + +++ + Unda, 


in which the p's are the respective probabilities of all possible different values 
of X. 


DEFINITION 2: The words “standard error of X”’ and the letter c, will be used 
to denote the expression 





oz = W/E[X — E(X)P. 


It will be noticed that the standard error of a random variable X may have 
its value equal to the standard deviation of some numbers u but that this does 
not mean that the two conceptions are identical or even similar. The E(X), 
and consequently oz, can be calculated only when the probability law of X is 
known, and are constant for the population from which samples are drawn. 
On the other hand, S can be calculated for any sample of the population and 
changes in value from one sample of w’s to another. 
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Before proceeding to the third conception, that of an estimate of the standard 
error, which is occasionally confused with the standard deviation or the standard 
error, the unbiased estimate of a parameter must be defined [3]. 

Consider a set of n random variables X,, X:,--- ,X,. These may be, for 
example, the volumes of timber to be observed on n strips that are going to be 
selected from some area by one random method or another. Denote by @ a 
parameter involved in the probability law of the X’s. For example, 6 may be 
the total volume of timber in the area. 

Let F be any function of the X’s. 

DEFINITION 3: If it happens that the mathematical expectation of F is 
identically equal to 6, then it will be said that F is an unbiased estimate of 0. 

Usually there will be an infinity of unbiased estimates of a parameter 0. 
They may be classified by the nature of the function F. Thus linear estimates 
may be considered such that 


F = Xo + MX + es + AXn 


in which the )’s stand for some fixed numbers. 

DEFINITION 4: It will be said that a linear unbiased estimate of @ is the best 
linear unbiased estimate (B. L. U. E.) if its standard error is smaller than or, 
at most, equal to that of any other linear unbiased estimate. 

It happens frequently that, while it is possible to determine the best linear 
unbiased estimate F of a parameter, it is not possible to calculate the value of 
its standard error, cy. For this purpose it would be necessary to know the 
whole population sampled. In such cases an unbiased estimate of the square, 
oy, is calculated. An unbiased estimate of the square of the standard error 
of F will be denoted by »>. This is the third of the conceptions mentioned 
above. 

The reason for the extensive use of the linear unbiased estimates and of their 
standard errors considered as measures of accuracy is the so-called Theorem of 
Liapounoff. Its content can roughly be explained as follows: If the variables 
Xi, X2,--+-,X~, are independent and the number n not too small, then the 
probability that F — 6 will exceed a fixed multiple of oc, is approximately equal 
to the probability as determined by the normal law. The above conclusion 
remains true whatever the probability distribution of the X’s that is likely to 
be met in practice and also in certain cases where the X’s are mutually dependent, 
for example, when they are determined by sampling a finite population without 
replacement [4]. 

The above conclusions do not apply to estimates that are biased in the sense 
of the above definition. . Also the standard error of such an estimate would not 
be a satisfactory measure of its accuracy. 


3. Description of data. Complete inventory data from the Blacks Mountain 
Experimental Forest, located in the Lassen National Forest, provide suitable 
material for testing the applicability of sampling theory to timber cruising. 











182 A. A. HASEL 






The timber is a virgin, all-aged stand, classed as pure pine type, with more 
than 90 per cent of the volume in ponderosa pine and Jeffrey pine. Most of 
the volume is in over-mature trees, i.e., trees over 300 years in age. The stand 
is considered to be fairly representative of the medium and the poor site qualities 
of the northeastern California plateau. 

With the exception of a few localities, all of the area was mapped as of uniform 
timber type according to the standards commonly used. Being fairly uniform 
also with respect to site quality, it may therefore be considered as a single 
stratum. Variability of stand volume from place to place within a stratum may 
be generally expected to be less, on the average than variability between places in 
different strata. Likewise, within a stratum, variability within compact sub- 
divisions may be expected to be less than average variability within the whole. 
Heterogeneity can therefore be controlled somewhat by subdividing the stratum 
into blocks and treating each block as a separate population. 

More frequently than not, in practice, volume estimates are needed both for 
the total timbered area and for separate working units or compartments within 
the area. In general, working unjt boundaries are defined by roads, ridge tops, 
drainage channels, and regular land subdivision lines. These working units 
can be taken conveniently as blocks, or if large enough, may be subdivided into 
two or more blocks. Such is the basis used for subdividing the area in the 
present study. 

The complete inventory data for these blocks are given in Table I. All the 
strips are 24 chains in width and extend in an east-west direction. The length, 
X, is given in 10-chain units, and the volume, Y, is given in units of 1,000 feet 
board measure. 


4. Method of estimation based on correlation between volume and strip 
length. The usual practice in sampling timber stand volume is to take measure- 
ments on plots or strips that are either regularly spaced or selected at random 
from all possible plots or strips within blocks. Oftener than not blocks are 
irregular in shape, and the number of plots along lines or the lengths of strips 
will vary. This variation introduces the matter of proper ‘‘weighting” in calcu- 
lating sample statistics. Such is the case in 15 of the 20 Blacks Mountain blocks. 

If we let Y; represent the volume on the 7th strip of length X; , with length 
expressed say in 10-chain units, and assume that the entire block contains a 
population of N strips, then the average volume to the unit of strip is B = 

N N 

» Y; 2» X;. Itis obvious that, if >> X; is known, and this is assumed to be 
true, the problem of estimating 6 is equivalent with that of estimating the total 
volume. The usual procedure of estimating is this: 

Out of the N strips within the block a sample of n is taken, giving n pairs of 
numbers selected out of the X’s and Y’s. Let us denote them by 


U1, Yr 5 Te, Y25°** 5 Ens Yn- 
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The ratio b = >> y; jf >. z;, is then considered an estimate of 8, so that the 
i=l i=l 


N 
estimate of the total volume in the block is b )> X;. 
t=] 

Our purpose now will be to study the above estimate b from the point of view 
of unbiasedness. In this paper it is assumed that the sampling of strips is 
purely random. To find out whether b is an unbiased estimate or not, its 
expectation must be calculated. This will be done in two steps. To begin 
with assume that the values 2 , x2, --- , 2, are chosen in one way or another and 
fixed. The value of b will then depend on the y’s only. It is possible that to a 
given value of x, say xz, , there will correspond just one value of y; in the block, 
but generally there will be several strips of the same length 2; with varying 
volumes of timber. The selection of any strip of this group to be included in 
the sample will keep the denominator of b constant, but will cause some variation 
in the numerator. The expectation of b calculated under the assumption that 
the z’s are fixed is 


(1) Eb | 21, 25-05 tm) = D Buel) / DO x, 


in which F(y; | 2;) denotes the expectation of y; calculated under the assump- 
tion that z; has a fixed value. Obviously E(y; | x;) will be what is called the 
regression function of y on xz, or of volume on the length of strip. 

It is safe to say that the graph of E(y| x) would almost always be rather 
irregular. On the other hand, it is known that the substitution of smooth curves 
representing the regressions for the true irregular polygons frequently gives 
results that are surprisingly accurate. Therefore it would not be unreasonable 


to use the assumption that E(y | x) can be represented by a polynomial of some 
moderate degree, 


E(y |x) = Ao + Ait + Age? +--+ + Az’. 


Substitution of this expression in (1) gives 


A D ti , Xs 
E(b| 21, 22, °*+ te) = +A + A +--+» +A, SH. 

2% D 2 > 2 

— i=l i=l 


But this is the conditional expectation of b, calculated under the assumption of 
fixed x’s, is only an intermediate stage in the calculations. We need an absolute 
expectation, calculated under the assumption that the z’s are selected at random. 
This gives 


> 23 v2 
(2) E(b) = AcE “m + Ai+ AE —— bs 64S 


D zi Dt aE? 


t=1 t=1 t=1 
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The value of 8 has the form of (2), except that instead of E(22x7'/Zz,) it contains 
~X7/=ZX;. Since in general the former does not necessarily equal the latter, 
for the unbiasedness of b it is necessary and sufficient that Ag = Ag = +--+ = 
A, = 0. This condition implies that the regression line of y on x is a straight 
line and passes through the origin of coordinates, 


(3) E(y |x) = Ais. 


Whether (3) is satisfied is a question of fact and can be authoritatively an- 
swered only by direct studies of regressions on some extensive inventory data. 
It may be noted also that in order to presume that (3) is usually satisfied, it 
should be established for a large number of areas. On the contrary, if a study 
of only a few areas shows that (3) is not true, then it would not be wise to take 
it for granted when attempting to make a sampling inventory of an un- 
familiar area. 

To investigate this point, linear regression equations of volume on the length 
of strip were calculated for 15 blocks of the Blacks Mountain Experimental 
Forest and it was found that the constant terms were both positive and negative 
with their absolute values varying from 12 to 677. The conclusion drawn is 
that the usual estimate b of the average volume per unit of strip is likely to be 
biased and that there is justification in looking for an alternative method leading 
to unbiased estimates. 


5. Best linear unbiased estimate of volume, based on the linear regression of 


volume on length of strip. In this section will be suggested a method of es- 
timating the total volume, say @, of a timber stand, which could be considered 
as an improvement on the one considered above. The new method consists of 
using a linear unbiased estimate of 6. In order to deduce the form of this 
estimate, certain assumptions have to be taken for granted concerning the 
timber stands to be sampled, and if it happens that these assumptions are 
unsatisfied in a particular case, the new estimate will not necessarily possess the 
desired property of unbiasedness. 

In deducing .the estimate F it will be assumed that the timber stand to be 
sampled satisfies the following conditions: (1) That the regression of timber 
volume on length of strip, X, be (approximately) linear and (2) that the vari- 
ability of the ¥’s for a fixed N is precisely known. It will not be assumed, 
however, that the linear regression line passes through the origin of coordinates, 
and this will allow F to be unbiased in such cases, as exhibited above, where b is 
biased. Following the Markoff method [5], [3] it can easily be shown that there 
is an infinity of linear estimates of 6 which are unbiased under condition (1). 
It follows that a choice can be made among them so as to diminish the standard 
error. This, however, is possible only when something is known about the 
variability of the Y’s when the value of X is fixed. For the present we shall 
assume condition (2) concerning this particular point, but in practice this will 
generally be quite impossible. This point will be considered further in Section 6. 
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Consider then a sure or non-random variable’ X able to assume the particular 
values X,, X2,---,X,. Assume that there is a finite population 7; of N; 
numbers ui, Ui2, °** , Win, COrresponding to each value X;, 7 = 1, 2,---,s. 
Assume that the mean u;. of the population 7z; is a linear function of X;, i 
for any 2, 


e., 


= A + BX;, 


with some unknown values of A and B. 

Assume that out of each population z; there is selected without replacement 
a random sample of n; individuals, with 0 < n; < N;, and denote by yi, 
Yi2,°** , Yin, the values of the u’s to be drawn. 

If the regression of the amount of timber on the length of the strip is linear, 


then the problem of estimating the total stand is equivalent to that of estimating 


6= > Ni(A + BX) =4 DN +BONX,. 
i=l 


t= 


Since the length of the strip, X, could be measured from any arbitrarily chosen 


origin, no generality will be lost by assuming > N;:X; = 0, so that @ = 


i=1 

A > N; = AN (say). Weighting the y; ; equally for each fixed 7 the B. L. U. E. 

t=1 
of 6 may be denoted by 
(4) F = Do mryjyi- 

i=l 

in which y;. =.Zy;;/n;. Here the \’s must satisfy the conditions of unbiased- 
ness, 
(5) E(F) = 6, 


and of optimum, 
2 — 
(6) or = minimum. 


It may easily be shown that condition (5) will be fulfilled by (4) if the \,’s are 
so selected that 


(7) b is nur = N; y niri Xi = 0 


t=l1 t=] 


8 This is an English translation from an excellent French term ‘‘nombre certain’”’ and 


“fonction certaine’’ to denote a non-random number and non-random function, invented 
by Fréchet. 








188 A. A. HASEL 








Condition (6) may now be considered. From the general formula for the 
variance of a linear function of several random variables and the fact that y;; 
is independent of yx: , 








= D> {nriSi + nln; — IAEl(ya — us.) (ya — wi ]} 


t=1 


‘ 2 
(8) => | man Si - ie 71 (nid; — rand | 


i=l 
= > 8 a me Mi = Do Aza; (say), 
i=l t=] 
in which S; stands for the (S.D.)’ of the population 7;, i.e., 
Nix 
S; = & 7 (ui; = u;.)*. 
N; j=1 


In addition to satisfying equations (7), the \,’s must be selected so as to mini- 
mize (8). 
Using the method of Lagrange, we find 
















(9) di = i (a + BX), 


for the case where 0 < n; < N, and A; # 0. If n; = N;, then A; = 0 and 
a+ BX; = 

Assume first that all n; < N;,7 = 1, 2,---,s. Then a and @are obtained 
from equations (7) after substituting in them (9), namely 


a ke wu; +8 ~ wAX; = N 
(10) 1 z 
ws i Ww; A; = B Zz w; X? => 0, 


where, for simplicity 


2 é 

Nn; (N; — 1)n; i 
11 , = = (tcl ; > , = W . 
= . Ae (Ni; - nS? =” 








If w; is considered as the weight of the observations at X = X;,, it will be con- 
venient to introduce a weighted mean and weighted S.D. of X’s as follows: 


, R Wy « i : W; x 
= _ tml i v2 tm =2 

(12) i= Ww i= "= a r. 

With this notation equations (10) can be rewritten and easily give 

Ni 
|? — ~ Ws? 
| N(S; + 2°) 
a= =" 6 
| WS? 





(13) 
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Substituting these values into (4), simple transformations give 
(14) F = N(g — Zbo), 
in which 7 = } wiy:./W, and bo represents the unbiased estimate of B and is 
given by 
bo = [(1/W) 2X wi Xiyi. — 29)/S2. 

The next step is to calculate o,. Substituting (9) in (8) and using (11) and 

(12) gives 
oy = Wla + Bz) + We’S?. 

Using (13) gives finally 


Nn’ z 


If X is the length of a given strip in any chosen units and X the average of 
such X’s for a given block, then (14) and (15) may be written 


| F = Nig + bX — 3)] 
r2 7 _ =\2 
Similarly for the case where one of the n,’s equals N; , for example, nm, = Ni, 
we find 
(17) F = N[y. + (X — X)I, 


in which 


(16) 








ie > wi(X; — Xi)(yi. — wad 
2 w(X; — X,)? 
Also 
N*>(X — X,)’ 


(18) o, = . ; 
2 w(X; — X,)’ 

It should be emphasized that X, in the above formulae does not necessarily 

represent the smallest of the X’s but the one of them for which n; = N;. 

The case where two or more of the n,’s are respectively equal to the corre- 
sponding N,’s need not be considered in detail. Together with the assumption 
of a strict linearity of regression such an assumption, for example, that n; = Ni, 
and n. = N2, would lead to the conclusion that the regression of volume on the 
length of strip is accurately known and that the estimation of @ could be made 
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without error. Owing to the fact that the hypothesis about the linearity of 
regression is, at best, only approximately correct, the errors of estimation will 
always be present and it is imperative either to arrange the sampling so as to 
have at most one of the n,’s equal to the corresponding N;, or to base the 
statistical treatment of the sample on a theory different from the one con- 
sidered here. 



























6. Additional hypotheses concerning S?. The formulae (16) with the w,’s 
determined by (11) are impossible to apply in practice because we do not know 
the values of the S;. The best we can do is to make plausible guesses as to 
what may be the values of the S?. These guesses are bound to be at most 
approximately correct and therefore the estimates of 6 that one can apply in 
practice will be only ‘“‘approximately best.’’ It is easy to see, however, that we 
may keep them unbiased. 

Suppose that we denote by 7? the presumed value of S;. Substituting this 
value in place of S; in (8) we should repeat all the calculations, leading us to 
such \; that will assure the unbiasedness of, say 


8 
/ 
F, = Do marjgyi., 
t=1 
and also a minimum value of, say 


2 <9 ni(Ni 
i = » rT; a 
The values of the ; will be obtained from the same formulae as those of ), , 
except that instead of S; they will depend on 7;. Consequently F, will have 
the same form as F, 
(19) F, = Nig’ + bo(X — 2’), 


with the difference that z’, 9’, S., and bo will now have to be calculated with 
different weights, say 


ni) 


elie 
ao a. 





_ (Ns = 1m .*% 
salons (N; — m)ri hes 2. 





If the form of the unbiased estimate F, is as that of F, the square of its standard 
error 9, is more complicated. In order to calculate it we have to go back to (8) 
and substitute into it the new values of \; obtained from the guessed weights v:. , 


— N;- 1 ‘ 'X,) 
(Ni - nari sins 


with 


N ’ , 
a’ — Vs? (Ss? + I *) 
(20) 
_ _ Nat 
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we have 


a? s? me _ = Vg 
(21) tl 


= dX v; pila’ + p’X,)’, 


where p; = Si/r;. It will now be helpful to introduce notation for another 
kind of weighted mean and weighted (S.D.) of the X’s, with weights equal to 
vip;. So let us write 


dX ™ piX ” » M6 piX " 
(22) Zz" = SS et tego. oo OE, 


Soe U5 Pi » U5 Pi 
‘ 


Expanding (21) and using (20), we have 


N’ >» Mi pi _— a\2 xg !2 Ql 
7 Z'(z’ — Z”) z°S,; 
a) da i eT + a}. 


Formula (23) refers to the case where the X’s are measured from their popula- 
tion, mean, X. In order to reduce it to the case where the X’s are given in 
their original values we have to substitute (z’ — X) for Z’ and (#” — X) for 2”. 
Thus 


N? Dd vip; Ud * 
ay ga ef ae HT ae 





Applying a similar procedure to the case where m = N; = 1, but n; < N; 
fori = 2, 3, --- , s, we easily find 


| al - v(X; — X1) (yi. 7 Yi.) 
(25) yn EN a, Cg RE renee by 
| 2d vi(X; iam x,’ 


and 


> V; pi (X; = x,)’ 
(26) = N*(X, — X)* 32 


[Eo mh} 


This formula will help us to test the appropriateness of guesses about the 
values of S?. It will be noticed that the \’s contain S? or r; in the same powers 
in the numerator and in the denominator. It follows that all we need to guess 
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is a system of numbers proportional to S? and not the S? themselves. Our 
problem will be to test a few such guesses on the data of the Blacks Mountain 
Experimental Forest and see which of them gives generally a smaller value of o; . 

Table II gives values of the S;, calculated for 15 blocks, together with the 
corresponding X;. In a few cases N; = 1 and consequently S; = 0. These 
cases are not included in the table. Using the values of S? from Table II and 
assuming systems of the n,’s, the values of «7 were calculated for these blocks. 
These would be the true (S.E.)’ of the best linear estimates of the total timber 
volume in each block, but it would never be possible to calculate them from 
sample data. 

The o;’s were calculated using the following guesses concerning the S? : 
(1) That they do not depend on X;, (2) that they are proportional to X;, 
and (3) that they are proportional to »/X,;. The ratios o+/o% for all blocks 
taken together were found to be .770 for guess (1), .769 for guess (2), and .777 
for guess (3). It is seen that, on the average, the guess that the S? are propor- 
tional to X; gives the smallest average value of o,. It is interesting, however, 
to note that the differences between the three guesses are, for all practical 
purposes, negligible. 

Ratios like 0/0; are sometimes described as the “amount of information” in 
F, as compared with that in the best linear unbiased estimate F. This ex- 
pression was introduced by R. A. Fisher. In certain cases it has the following 
property which justifies the term used: Let n be the size of the sample which 
serves for calculating F, , then, if it were possible to calculate the best linear 
unbiased estimate F, the same accuracy of estimation would be obtained by 
using a smaller sample size no¢/o,. In the case considered in the present paper 
the above circumstance does not occur. Still, the ratio o;/o; seems to be con- 
venient to describe the situation. 


7. Another scheme for estimating 9. It will be noticed that the ignorance 
of what are the Sj is not the only circumstance which makes it difficult to apply 
the above formulae. There is also another one connected with the values of N; . 
We have N; = 1 in several blocks and for several different strip lengths. True 
this might have been avoided by defining block boundaries in such a way that 
N; > 2, but it was considered best to conform strictly to the practical situation 
where the N,’s may be smaller. In such cases we may include in our sample 
all the strips of a given length, say X,. If we apply to such samples the above 
formulae, deduced under the explicit assumption that the regression of Y on X 
is strictly linear, we shall force the fitted regression line through the point 
(X,, Y:.). As the assumption of strict linearity is obviously not exact and the 
exhaustion of strips of length X, is possible only when there are very few such 
strips, the whole procedure may lead to serious inaccuracies in the final estimate. 
One safeguard against this is never to exhaust strips of any given length when 
dealing with formulae deduced from finite populations. 

The fact that the true regression point (X,, Y,.) does not actually lie on a 
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straight line makes it uncertain whether taking into account the finiteness of 
populations of strips of the same length is beneficial to the accuracy of the 
finite estimate. In the preceding sections we worked on the assumption that 
there is but a finite number of strips of the same length and on an inaccurate 
assumption that the regression is strictly linear. In the present section the first 
assumption will be dropped, having in mind that the effect of the inaccuracy of 
the second assumption may thereby be reduced. 

The assumption that each of the N; is infinite will be made only in deducing 
the A; and will be reflected in weights. Formula (11) will now reduce to w; = 
n;/S?. If we assume further that S? = X7/k, where y and k are some con- 
stants, then 


kn; nN; 
= > W=Lm=kd x, 





t t 












and the final estimate is 


(27) F = Nig + 6(X — 3). 
The square of the standard error of F has again the same form as in (16), 
2 _N’ (X — 2)’ 
- on ae "ss fF 


the only differences being in W, , and 82. If y = 0,so that the S; are assumed 
to be constant, then 


DO; = kn; ; W=kdin., 


and all the symbols z, 7, and S? assume their customary meaning of ordinary 
means and of ordinary (S.D.)’. 

It would be easy to deduce explicit formulae for y = 1/2, etc., but they are 
not elegant and, if the necessity arises, the calculations could be carried through 
by starting with ®; = 1/X]7. The omission of k does not influence the form of F. 

The question whether the combination of one true hypothesis about the N; 
being finite, with another incorrect one that the regression is strictly linear, is 
better than that of two incorrect hypotheses, will be studied by means of a 
sampling experiment in Section 9. 


8. Unbiased estimates of c;. While it may not be unreasonable to hope 
that a guess of a system of numbers proportional to the S; may be successful, 
it is entirely hopeless to try to guess the actual values of the S;. It follows 
that, if it is desired to obtain from the sample some sort of measure of the 
accuracy of F, we have to calculate an estimate of o7 . 

We shall treat the problem by assuming that the regression of Y on X is 
strictly linear and that the Sj are proportional to X7 and that the N; are all 
finite. It will be noticed that they will enter the formulae by means of the 
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ratios (V; — 1)/(N; — n,). If it is desired to obtain formulae referring to the 
assumption of infinite N,’s, it will be sufficient to replace these ratios by unity. 
Of course, the symbol N will always represent the total number of strips in the 
actual block on which it is desired to estimate the volume of timber and will 
not be affected by the assumption of the N,,’s being infinite. 

On these assumptions 


E(y;.) = A + BX,, 
Sue = E(yi;— A— Ow me St = kX}, 
with some value of y supposed to be accurately guessed, which however need 


not be specified, and with an unknown factor of proportionality, k. Thesquare 
of the standard error of y;. is then known to be 


f _* i) 
(29) é.. ae S; N; nN sa XI(N; ) 


Ns N; —1 n(N; — 1) p 





The right-hand member of (29) is equal to the reciprocal of what we have 
formerly denoted by w; and described as the weight of the observations at 
X = X;. We have mentioned above that the formula giving F does not 
depend on the values of the w; , but on proportions between the w;. In other 
words, if we drop the unknown factor k and denote by w; the ratio 


X7(N; — ns) 


which involves only known quantities, these new weights will lead to exactly 
the same value of F as the original weights. It will now be convenient to alter 
the definition of weight and use formula (30). With this new meaning of w;, , 
(29) could be rewritten o),. = k/w;. 

Let us further use the letter m to denote the number of those X,’s for which 
we have at least one observation. In other words m will be the number of 
different lengths of strips in the sample and also the number of different y;.’s 
that we are going to calculate from it. 

Now let us go back to formula (16) giving the square of the standard error 
of F. We notice that, while and ‘S? do not depend on the unknown factor of 
proportionality, k, the sum W of the original weights does depend on it and with 
our new meaning of w; , 





1 
W = k » Wi. 
It follows that o> should now be written in the form 
kN’ l (xX — =) 
31 - 144-77) 
( ) Or z w; + Ss? 


s 
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and that, in order to estimate a7 it is sufficient to get an estimate of k. We 
easily get an unbiased estimate of k by merely applying the second part of the 
Markoff Theorem [3]. According to it an unbiased estimate of k, based on 
m — 2 degrees of freedom is given by the ratio 


(32) oe 9 — Ne — OF 


t 
i=] m— 2 , 


in which Z, 7, and by are calculated according to the assumptions made regarding 
N,; and y. It may be expected, however, that the estimate (32) will not be a 
very accurate one because the number of degrees of freedom on which it is based 
may be very small. 

In an attempt to find a better estimate of k we shall proceed by analogy and 
calculate the expectation of a sum similar to the one in the numerator of (32) 
but depending explicitly on the particular y;,’s, namely of 












S= >> Ws — 9 — bX, -— OP. 
i=l jel ny 
It will be noticed that if the N; are finite, y;; and y;; are dependent and that the 
Theorem of Markoff does not apply to S$. Introduce the notation 


ni = Yi — A — BXi, 


1< 
nN = —- > ni = y;— A— BX;. 


Ny j=1 
Easy, but somewhat long calculations show that S$ can be rewritten in the form 


m 2 1 m 2 
—" (= win) + 3 P w(Xi — Dn 
Ss? = Zz 7. Ws na — \i=1 z Li=l , 


i=] j=l Ns — 
Wi 


t=] 










which is most convenient for calculating the expectation sought. We notice 
first that 


E(nii) = kX?, 





2 2 c 
E(n;.-) ee. — =< 
Wi 


Further, as y;. and y;. are mutually independent if 7 # j, the same is true for 
ni. and n;.. It follows that 


E(yi.n;.) = 0, tj. 


Consequently 


m 2 m m m 
E(S wn.) = BZ wtet.) = Det) = kd ww. 


i=l] i=] t=1 
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Similarly and for the same reason 
E(X w(x; — Z)n:.]? = ks? ie Wi. 
It follows that 


E(S?) = e| 2 AD a 2|, 


1 N; — 7; 
and that the ratio 


(33) So DY low — 9 — boxe — Pe 
33 <gremeemmrentonmaeeamennes GR Seabsetcagrnmmenmpmencencenmemmnannentent: 
~~ = ns — 1) 
EM m 7 Nam? 


is an unbiased estimate of k. In cases where either all n; = 1 or all N; are 
infinite the denominator of (33) reduces to the number of degrees of freedom 
in So, equal to =n; — 2. In other cases the denominator of (33) is greater 
than the number of degrees of freedom. Whether the numerical difference is 
large or small depends on the fractions (V; — 1)/(N; — n,;). We may expect 
that in many practical cases it will be small. 

We shall write 


| 


S, = YZ wo / Dm, 


tml 


t 
\ 
t 
e 
ri 
f 
tz 
‘e. 
S 


dX w(X; — Z)y:. 2d Ww; 
S. Sy ; 


T= 


ee = Sees 
* vere S 


It follows that 


Si = : wSi(1 — 7’). 


etree 


Substituting this formula into (33) and then the result of this substitution in 
place of k in (31), we finally get 


2 2 =\2 
(34) ee ee fe E +G-2 |. 
_ n(Ni— 1) _ 2 S: 
i= Ni-n; 
The case where one of the n; is equal to the corresponding WN; , e.g., where 
n, = N, = 1 is treated in a similar manner. Using formula (18) and the nota- 
tion adopted above, we can write 


N*(X, — X)? 


2 

Or = k - 
Z. wi(X; = X,)’ 
ton 
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The unbiased estimate yu; of o will differ from this expression in that instead 
of the unknown factor k it will contain its unbiased estimate. To find this 
estimate we proceed exactly as above and calculate the expectation of 


So = 7 7 [yi = je — bo X; = Xy)/ =. 
im? jul Ns 





2d wi(X; -_ X1) (yi. = Yi.) 
bo = a ea i Ei a ei 


a wi(X; a x,)’ 
The unbiased estimate of a7 is 
(35) 2 Ss N?(Xi = xX)’ 


fe - >= — : 
Po Ws F wk. ~ 2S 


imz N; — 1; i=2 


The number of degrees of freedom, f, on which yu? is based is 


(ian mt, 


t=2 





9. Empirical tests of the preceding theory. Applications of any mathematical 
theory involve certain assumptions about the phenomena studied that are not 
exactly true. In order to have a reasonable hope that the predictions of the 
theory will be comparable to the actual facts, we must perform empirical tests 
and see whether such deviations from the assumptions underlying the theory as 
are usually met in practice influence materially or not the working of a given 
theory. Our object in the present section will be to test whether and to what 
extent such deviations influence the applicability of the theory. For that pur- 
pose it will be useful to enumerate the more important uses of the theory that 
are likely to be made. 

The first point refers to the choice of the standard error o, of the best linear 
estimate F as the measure of accuracy with which F estimates the unknown 
volume of timber, 6. If all the assumptions were true, the Theorem of Lia- 
pounoff would guarantee that, when the size of the sample, =n; , is only mod- 
erately large, the frequency distribution of the ratio 


(36) (F — 6)/or, 


would be very approximately normal about zero with unit S.E. If this were 
actually true then the value of or would be a justifiable basis for the choice 
between various alternative estimates of 6. However, the discrepancies between 
the hypothesis underlying the theory and the actual facts may easily produce a 
bias in F, or may deprive or of the above important property. 
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Therefore, the first thing that we have to test is whether in such conditions 
as are actually met in practice the ratio (36) is in fact distributed in repeated 
sampling in a way that is comparable with the normal law. The data of the 
100 per cent survey of the Blacks Mountain Experimental Forest will serve us 
for the test. 

The second important application of the theory is connected with the use 
of uy. The purpose of calculating yu, is to characterize the accuracy of the value 
of F obtained from the sample. The most appropriate way of doing so is to 
calculate the confidence interval for 6. This has the form [5] 


(37) F — tour SOS F + tap 


in which ¢, denotes the ‘‘Student’’-Fisher ¢ taken in accordance with the number 
of degrees of freedom in yu, and the chosen value of P. The confidence interval 
has the property that, if calculated for a great number of samples, the frequency 
with which the true value of @ will lie between the limits F + t.ur will approach 
the value a = 1 — P defined as the confidence coefficient. 

The above statement concerning the confidence coefficient is strictly true if, 
apart from the various hypotheses that were enumerated, the distribution of the 
y’sisnormal. As a result of a theorem by Kozakiewicz [6] the same statement 
will be approximately true also for non-normally distributed y;,’s, on condition 
that the sample sizes are considerable. In the situation where the above theory 
is to be applied all these assumptions are not satisfied. Still the formula for the 
confidence interval may well work, but before accepting this we have to have 
some experimental evidence. The crucial point that it must cover is whether 
the ratio, say 


(38) t = (F — @)/ur, 


does or does not follow in repeated sampling a distribution which is sufficiently 
close to the theoretical one, known as “‘Student’s” distribution. If the empirical 
distribution of t does approach ‘‘Student’s”’ law, then the frequency of correct 
statements concerning @ in the form (37) will be approximately equal to the 
chosen a, and conversely. 

The n,’s for this experiment were fixed according to the systems shown in 
Table III, with all X’s having a chance of appearing in the samples, and the 
n,’s quite closely proportional to the N,,’s and approximately 25 per cent of the 
latter. Random sampling numbers [7] were used in making the selections of n; 
strips out of any group of strips. A total of 150 block samples were drawn, 
equally distributed among the 15 blocks. 

There are 95 samples for the case where alln; < N;. For these, formula (19) 
was used to calculate F and formula (24) for o+ , using the guess that the S; 
are constant over all strip lengths. On the hypothesis that the ratio (36) is 
normally distributed about zero with unit S.E., we divide the range of variation 
of possible values of (36) into 20 intervals such that, if the hypothesis is true, 
then the probability of an observed value falling in any particular interval is 
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TABLE III 


Systems of n;’s for sampling experiment 


Block 


9 











Total 





Xi 














Total 


20 


Total 


equal to .05. For 95 samples then, the expected frequency in each interval is 
4.75. The observed frequencies are shown in Table IV. 
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The agreement between the observed and the hypothetical distribution is 
tested by means of the fourth order smooth test for goodness of fit [8]. The test 
is designed so as to be particularly sensitive to such deviations from the hypo- 
thetical distribution that could be described as “smooth.” It is used here 
because it is expected that, if the actual distribution of the ratios considered 


TABLE IV 
Frequency distribution‘ of (F — 0)/or and r= 6)/ur calculated under various assumptions 


Assumption of finite population of strips 





| Assumption of infinite population of strips 


















































| 


= (F — O)/ep (F — 0)/up (F — 6)/up 
All nj < Ni One j= Ni | All ms <M | one m= Ny || All ni < Ni =? ‘ Total 
ke we a nk | nk | om nk 

5 11 4 9 | 3 4 7 

3 1 5 1 4 | 2 6 

5 2 6 2 7 2 9 

8 1 4 0 2 } 4 3 

3 2 4 2 |} 5 | 4 9 

4 2 5 3 | 4 3 7 

8 1 6 2 | 6 | 0 6 

3 2 | 4 2 | § | 4 9 

3 2 5 2 | 8  ¥ | 10 

7 0 6 1 | 5 | 3 | 8 

1 0 | 3 | 0 | 4 | 0 4 

4 1 | 3 0 | 4 4 8 

6 1 5 2 7 5 12 

5 1 6 3 3 | 7 | 10 

5 1 6 1 | 7 _ rr 

_ i — — 

2 1 6 3 | 9 | 5 14 

5 2 8 4 | 6 ; sf 7 

4 6 5 | 3 df | 0 4 

10 1 3 | 2 2 | 6 . 

4 6 1 2 0 } 1 1 

Total 95 44 95 44 | 95 | 55 150 

vi: 1.326 33.812 | 5.463 | 13.091 | 9.055 2.572 | 8.764 

P 85 <.0l | 25 | O1 || 06 63 | 07 
P(x?) .57 <.01 | 63 09 || 18 12 | 21 





does differ markedly from the normal or from “Student’s’” one, then still the 
curve representing this actual distribution would be a “smooth” one, presumably 
with a single mode, and would cross the hypothetical curves at only a few 
points. There is empirical evidence to show [9] that in such cases the smooth 
test of fourth order is more powerful than the usual x’ test. 


4 By 20 intervals of equal probability. 
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The criterion used in the smooth test of the fourth order is denoted by yj. 
If the hypothesis tested is true, then yj is distributed, approximately, as x’ 
with 4 degrees of freedom. To calculate yi we proceed as follows: 

Let x be a random variable and H denote the hypothesis that the distribution 


of x is given by a perfectly specified function f(x). The range of variation of x 
is divided into 2s = 20 intervals, 


(— oO, a), (a, , a), ne (ayo , +), 


so that, if H is true then the probability of x falling within any such interval is 
exactly equal to .05. Such a subdivision can frequently be made easily from 
appropriate tables for f(z). We associate with these intervals a variate z whose 
value corresponding to the kth interval will be 

_2k-1 1 _ &Ak—s)-1 


oe? ota Fh ww 8 oo 
” 48 2 4s , sini ase 


It will be seen that if we start at the point a, and follow up the intervals to the 
right and to the left, then the corresponding values of z will be 

—s 3 a~ i 
2 — +7’ a 7 + 4s 7 
Consideration of the variable z is then substituted for that of the observed 
values 2; , %2,°°-,2, Of x. If any value z,, falls in the kth interval a, < 
Im < a, then this is interpreted as an observation of x which yielded the 
value z,. Let n, denote the number of observed x’s which fall in the interval 








28 
(ax-1, a.) and let the Gaussian symbol [z‘] stand for the sum [z‘] = >> mez. 
kewl 


To apply the fourth order smooth test such sums have to be calculated for 
a = 1, 2,3, 4. Then they are substituted into the equations below, deduced 
under the assumption that the number of intervals of subdivision of the range 
of z is equal to 2s = 20. 


uy = n-’(3.468,440[z])”, 

n *(13.500,884{z2"] — 1.122,261n)’, 
n'(53.857,548[2"] — 8.031,507[z])’, 

uy = n-'(218.148,007[z‘] — 46.239,587[z"] + 1.139,500n)*. 


Finally yf = uj + ui + uj + ui. If the calculated value of yi exceeds the 
tabled value of x’ with four degrees of freedom, corresponding to the chosen 
level of significance, then the hypothesis tested, H, should be rejected. 


2 
Ue 


2 
Us 


5 The above expressions for the u’s differ a little from those published in the original 
paper on the smooth test because in the latter the test was designed to apply only to un- 
grouped observations. The present formulae obtained in the Statistical Laboratory of the 
University of California appear in print for the first time. Obviously if the number of 
intervals 2s is increased, the formulae for grouped data will approach those for un- 
grouped ones. 
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The agreement between the observed distribution and the expected distribu- 
tion is shown to be excellent in Table IV, the probability of a greater difference 
occurring through errors of random sampling alone being .85. The correspond- 
ing P for the x’ test, where consecutive pairs of intervals are combined to make 
10 intervals in all, is .57. 

For the case where one n; = N; = 1 there are 44 samples. For these samples 
the value of F was calculated from formula (25), and the value of o> from 
formula (26), again taking the values of the S? as being constant over strip 
length within blocks. In this case the deviation from expectation shown in 
Table IV is greater than can be attributed to chance alone. These results are 
also obtained by the x’ test, which gives x’ = 25.091 and P < .01 on 9 degrees 
of freedom. 

The conclusions we draw from these results where one of the assumptions 
made is that the population of strips is finite, are that the block boundaries 
should be so defined that all N; > 1, or if this is not done, that the systems of 
n;’s be such that no sampling is done from strips where N; = 1. The fact that 
some n; = 0 when the corresponding N; > 1 has no appreciable effect on the 
working of the theory. In the previously described test for samples in which 
n; < N; the N used in formulae (19) and (24) always referred to all strips in the 
block, regardless of the fact that strips of some specified lengths X; did not 
appear in particular samples. 

Using the same samples, the distribution of (F — 6)/ur is compared in a 
manner parallel to that described above, to the distribution of the ‘‘Student’’- 
Fisher ¢, taking into account the number of degrees of freedom. 

The formulae used for estimating o; , namely for calculating yu}, are (34) 
where all n; < N;, and (35) where one n; = N;. The estimates of 6, namely F, 
remain unchanged from those previously calculated. 

The results from this second application of the theory as judged by the smooth 
test in Table IV lead to the same conclusions as were made from the first applica- 
tion of the theory, namely, that under the assumption that the population of 
strips is finite no N; should be exhausted in the sampling. 

It is interesting to note that the application of the x’ test to the observed 
distribution of (F — 6)/ur corresponding to samples with one n; = N; = 1, 
did not reject the hypothesis that it follows ‘‘Student’s” law. In this case the 
range of ¢ was divided into 10 intervals of equal probability and the value of x’ 
obtained was 15.091. With 9 degrees of freedom this gives P of the order of .09. 

The ratio (36) cannot be determined under the assumption that the population 
of strips is infinite where one n; = N; because the values of S?/r; cannot be 
obtained for such strips. Under this assumption it is impossible to calculate 
or by the formulae deduced in the present study and the first use of the theory 
must be omitted. However, the estimate of o} from samples can be calculated 
and the ratio (38) determined. 

The estimates F were calculated using formula (27), taking n; = w;. This 
same formula applies whether or not one or more of the N; are exhausted. Each 
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sample from Block 15 and one sample from Block 12 exhausted two or more 
strip lengths and their estimates could not be calculated under the assumptions 
made heretofore, but these can now be obtained under the present assumptions. 
The estimates u} were obtained from (34) for all samples, taking the S? as 
constant over all strip lengths and n; = w;. The fact that one or more N; 
are exhausted does not change the procedure for such samples in any way. 

For the case where all n; < N; in Table IV, the value of P = .06 obtained by 
the yi test indicates that the agreement of the observed distribution with 
expectation, although not close, is acceptable. When the data are regrouped 
into 10 classes and the x’ test is applied, we get P = .18 on 9 degrees of freedom. 

The yj test applied to the distribution of (F — 0)/u» for samples where one 
or more n; = N; indicates that the correspondence with expectation is good. 
This result is in marked contrast to the corresponding results in previous tables 
and bears out the belief previously expressed in Section 7, based on intuitive 
considerations, that by dropping the assumption of finiteness of number of 
strips of a given length, the error of the assumption of strict linearity of regression 
would be compensated for to some extent. On the basis of these findings we 
can add the conclusion that if, in sampling, the number of strips of a given 
length are exhausted, the assumption of finiteness should be dropped and the 
sample estimates calculated from formulae deduced under the assumption that 
all N; are infinite. 

There remains some question as to statistical treatment of samples in which 
alln; < N; , that is, whether to use formulae deduced for finite or infinite popula- 
tions. The final choice can best be based on the relative size of the confidence 
interval (37). Where all n; = 1 the estimates are the same under both as- 
sumptions. For estimates of all blocks taken together the finite population 
estimates tended to be within 5.5 percent of @ in 95 out of 100 trials, while the 
corresponding percentage for infinite population estimates was 6.0. We there- 
fore conclude that it is better to use the assumption of finiteness of N; where all 
ni < N;. 

The method of sampling considered here is what could be called restricted 
random. The restriction consists in that we group together the sampling 
units of the same size, select nonrandomly several such groups, and only then 
proceed to draw at random n; units of a group of N;. Frequently the strips 
of the same size will be situated within the block close to one another. In those 
cases the restricted sampling considered will assure that the sample will contain 
elements more or less uniformly distributed over the area of the block. 


10. Summary. Several methods of sampling timber stands and statistical 
treatment of the samples were considered. Data from a complete inventory of 
the Blacks Mountain Experimental Forest served for testing the methods in 
practice. 

It was found that the usual method of estimating from strip samples taken 
within nonrectangular blocks of timber gave biased estimates, unless the linear 
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regression of volume on strip length passed through the origin of coordinates. 
It was shown that this condition was not a safe one to assume. Consequently 
methods of estimation were sought which were freed from this restriction. 

The appropriate formulae for the best linear unbiased estimates were deduced 
under various combinations of the following assumptions: 

(1) That the regression of timber volume on strip length is strictly linear, but 

may or may not pass through the origin of coordinates. 

(2) That the values of the (S.D.)’ of timber volumes on strips of equal lengths 
are (a) constant for different strip lengths, (b) proportional to strip 
length, and (c) proportional to the square root of strip length. 

(3) That the number of strips of a given length in each block is (a) finite, 
and (b) infinite. Assumption (b) was based on intuitive considerations 
which indicated that this assumption, though known to be false, might 
compensate for another false assumption, namely, that of strict linearity 
of regression. : 

It was empirically found that assumption (b) of (2) gave better results than 
either (a) or (c). However, the advantage was small and, in the author’s 
opinion, did not justify the extra labor in calculations which are simpler when 
assumption (a) is made. Therefore all other calculations were made on that 
assumption. 

An extensive sampling experiment was made to test whether the smallness 
of the samples combined with the conflicts between assumptions of the theory 
and the actual facts, influenced the validity of the normal theory. 

Whenever the sample did not exhaust strips of a given length, it was found 
that the formulae based on the assumptions that the populations of such strips 
are finite and that they are infinite both work satisfactorily, generating distribu- 
tions similar to those determined by the normal theory. However, the confidence 
intervals based on the true assumption that the populations of strips of equal 
length are finite, proved to be narrower. Consequently, whenever the sample 
does not exhaust all strips of any given length in the block, the true hypothesis 
concerning the number of such strips should be used. Formulae (19) and (34) 
are therefore the appropriate ones, using weights based on finite populations. 

In cases where the sample did exhaust the strips of a given length, the treat- 
ment of the number of such strips as finite, combined with the inaccuracy of the 
assumption that the regression of timber volume on length of strip is linear, 
resulted in marked disagreement between the actual distributions of statistics 
and those based on normal theory. This disagreement was not found to exist 
in statistics calculated with formulae (27) and (34) used on the assumption 
of an infinity of strips of a given length. This suggests the conclusion that the 
exhaustion of strips of a given length by the sample should be avoided and, 
when this is impossible, then the formulae based on the assumption of an infinity 
of strips of a given length should be used. 

The formulae deduced can be applied equally well to line plots as to strips. 
With the formulae deduced the most efficient sampling will be obtained when 
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the average sample strip length is close to the average for the population, but 
where the variation among sample strip lengths is the maximum. 

The author is deeply indebted to Prof. J. Neyman for guidance and advice, 
and to Miss Evelyn Fix and Miss Phyllis Burleson for help in computa- 


tions. 
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TABULATION OF THE PROBABILITIES FOR THE RATIO OF THE 
MEAN SQUARE SUCCESSIVE DIFFERENCE 
TO THE VARIANCE 


By B. I. Harr 
Ballistic Research Laboratory, Aberdeen Proving Ground 


with a note 


By JOHN VON NEUMANN 


In recent publications von Neumann has determined the distribution of 
5'/s’, the ratio of the mean square successive difference to the variance, for odd 
values of the sample size n' and for even values of n.’_ In this paper the prob- 
ability function, i.e., the integral of the distribution, is evaluated for specific 
values of n. 

Let x be a stochastic variable normally distributed with mean ¢ and the stand- 
ard deviation o. -The following customary definitions for the sample are: 


._ 1 
the mean, se = Ze Ty» 
n p=1 
2 1 - 2 
the variance, c= = he (x, — =)’, 
p= 


2 


al 
and the mean square successive difference, 5 = — > (t.41 — 2,)°. Letting 
—-i st 


Fan 


s? n-—l1 





(1 — e), von Neumann shows that the distribution of ¢€, w(e), is 


symmetrical with zero mean and intercepts equal to + cos = (loc. cit.’, p. 372), 
and that w(e) is determined for odd values of n by 


a ye adie N-mr 


(n—1)—1 - _ 
dé™ + n ( =) 

I I € — cos — 

n 


in the odd intervals 





T 2n 
cos —- 2 e 2 cos—, 
n n 














cos 34 4 cos (n — 2)x cos (n — 1)x 
tekeet,-...2e-& ete ( a 

nm n n n 
1 John von Neumann, ‘‘Distribution of the ratio of the mean square successive differ- 


ence to the variance,’’ Annals of Math. Stat., Vol. 12 (1941), pp. 367-395. 

2 John von Neumann, ‘‘A further remark on the distribution of the ratio of the mean 
square successive difference to the variance,’’ Annals of Math. Stat., Vol. 13 (1942), pp. 86- 
88. 
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}(n—1)—1 
and by gana“) = 0 in the even intervals 
- 











an eekewn. 
et ncnent,--- nto scrc0mt—™ 
n n n 
(loc. cit.’ pp. 389-390). 
For n = 3, 
1 1 
for cos ~ = — 
or CO g =<s 3° 
For n = 5, 
") ao . 1 
GG) Saas ss: SSS 
w/—é+ie—ds 
(2) cont 4+-en0 pul ane 
2 a 5 5 5 s  § 
of, ae a, ss 
T T Tv 
6— a panel = oe ae 
cos + cos , COs = cos = COS = € cos 5 + cos — 
T 2r 3r Ar 
~ om oe = =. 
for cos= 2 €2 5 and cos 5 2e2 cos 5 
But for cos = Z2eZze Fw) = 0, thus 


(3) w(e) = const. 


For n = 7, 














(4) Manat ; 


rVf—&+5A—324+4, 








for cos 2 ¢e 2 cos = and cos => =e 2 cos - with the + sign, and for 
Z 
3a dr 
cos 2 e2 cos 7 with the — sign. 
But for cos a 2e2 cos 7 and cos = =e2cos *, w’’(e) = 0, thus 





(5) 


? The square of the modulus. The numerical evaluation of the inverse sine amplitude 
function used for n = 4, 5, 6, is taken from unpublished tables of the Legendrian elliptic 
integrals by F. V. Reno of the Ballistic Research Laboratory, Aberdeen Proving Ground. 
The square of the modulus is the argument for this tabulation. 


w’(e€) = const. 
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For even values of n von Neumann shows that the distribution of e, 


_ T3(n—1)] * €\ jn-3/, 4 “2 
wa+0) (€) — Tii(n — 21°) te Wa @) p (1 p) dp (loc. cit. ). 


cos — 
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For n = 4, 
Gs 
(e) = ——== 
mens vie PV1— p’ 
where w,{-—}) =— | — {-— cos-— }i -— cos— > 
T 4 4 


(6) dp 


as is 1 ee 
wa+ oy (€) Var Il, Vp — V2 ep + V26)(1 — p) 
oe — V3 wt (1,3- V3" 
wea eae) 
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forcos» >e>—~. 
’ inline 
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F = 6, = we 
or n wa+(o (€) i nel vi Vi =e 


oa sen a (oseolh ese] 2-64) 


dp, where 








for cos = = ¢ & cos * and cos &<¢¢ cos °™ , and where 


(8) w(e) = const. 


. 24 
for cos 5 zea cos >. 
The integrals needed to obtain w(¢) for n = 6 and w’(e) for n = 7 have been 
evaluated by numerical quadrature. Graphs of the distribution of 6°/s’, w(8"/s’), 
for n = 3, 4, 5, 6, 7, are shown in Fig. 1. 


The probability function, P(s’/s? < k) = [ w(3"/s*). d(6”/s*) has been ob- 


tained from w(8°/s’) by numerical quadrature for n = 4, 5, 6, 7. The results 
are given in Table III. 

As is mentioned by von Neumann, R. H. Kent has suggested a series ap- 
proximation of the form 


w(e) = } as (cos* © — 7; , 





since the order of vanishing of w(e) is 3n — 2, and since w(e) is an even function 
of ¢ (loc. cit.’ p. 391). Determining the a, by the condition of normalization 
and by the first three even moments of the actual distribution, M2, M, and Mg 
(given on pp. 377-378, loc. cit.’), and integrating the result, we obtain 
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tS eo 


9 


Ccos* 


L.(4[n + 21, in + - + 


T vi 
3 cos! — 15 cos® — 
n n 


_ MaBn + 190 +3) 4 —" ms 


4 (n+ Sn + (nm + 9) 1,(3[n + 4], 3[n + 4]) 5 _ Mn + 3) 


2 7 
oe — 
n 





4 M,(n + 3)(n + 5) _ Moe(n + 3)(n + 5)(n + "| 


Tv Tv 
3 cost — 45 cos® - 
n n 


The Tables of the Incomplete Beta-Function’ can be used to evaluate (9), with 

t= : (4 — + 1), Table I shows the results obtained for the eighth and 
2 \cos (2/n) 

tenth moments for the distribution (9) and for the true distribution for certain 

values of n. 


2 
Table II gives a tabulation of P(% < r) for n = 7 by the use of (9) and by 


the method of (4) and (5). The approximation (9) has been used for the com- 
putation of the probabilities of Table III for n = 8.° 

It has been shown (loc. cit.’ pp. 378-379) that for n — « the distribution of 
« becomes asymptotically normal. For n = 60 values of °/s’ are given below 
for different levels of significance. These values have been computed from 
Table III and from a table of the integral of the normal function with standard 


2n oy 
deviation equal to sg / aa iet 4)’ the square root of the second 
moment of the distribution of 8°/s’. 


4 Karl Pearson (Editor), Tables of the Incomplete Beta-Function, London: Biometrika 
Office, 1934. 

’ The results obtained by L. C. Young using the Pearson Type II distribution are suffi- 
ciently precise for the significance levels and sample sizes tabulated. Cf. L. C. Young, 
“On randomness in ordered sequences,’’ Annals of Math. Stat., Vol. 12 (1941), pp. 293-300. 
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TABLE I 
‘i Ms Ms Mw Mw 
(9) True (9) True 
é .00412 .00413 .00201 -00202 
8 .00318 .00318 .00150 -00151 
9 .00246 .00246 00111 .00112 
TABLE II 
8] 
P = <k)forn =7 
8°) 
k By (9) By (4) and (5) 
.25 .00001 .00001 
.30 .00007 .00007 
.o0 .00027 .00027 
.40 .00065 .00065 
.45 .00124 .00126 
.50 .00209 .00214 
.55 .00326 .00333 
.60 .00478 .00486 
.65 .00671 .00678 
.70 00911 .00913 
0 .01203 . .01197 
.80 .01552 .01534 
.85 .01964 .01932 
.90 .02443 .02403 
.95 .02995 .02957 
1.00 .03624 .03598 
1.05 .04333 .04325 
1.10 .05126 .05137 
1.15 .06006 .06036 
1.20 .06976 .07020 
Values of 8°/s’ for Different Levels of Significance 
n = 60 
P = Ol P= 005 P= Ol P = 05 
TE TEE. cc ccccccess 1.2558 1.3779 1.4384 1.6082 
eee ak 1.2358 1.3688 1.4333 1.6092 


This work was undertaken at the suggestion of Mr. R. H. Kent. I am 
much indebted to him and to Professor John von Neumann for many important 
suggestions and criticisms. 


Note to Fig. 1, by John von Neumann. Inspection of the graphs of w(8’/s’) 
for n = 3, 4, 5, 6, 7 (see Fig. 1) discloses certain singularities of the function 
w(é'/s’), which seem to deserve attention. 








TABLE III 


p(E<e)= [ o(2)e(2) 



























































; i | 4 5 | 6 | >i #} wo | wn 2 
| Pee per aa ee 
25 | | .00001 | .00001 | .00001  .00001 | 
30 | 00007 | .00007 ; .00005 | .00004 | .00002 | .00001 
35 | 00006 | .00027 | .00021 | 00014 ; .00009 | .00005 | .00003 
40 .00047 | .00065 | .00047 | .00031 , .00019 | .00012 | .00007 
45 | 00126 | .00126 | .00088 | .00059 | .00038 | .00025 | .00016 
.50 00038 | .00246 | .00214 | .00150 | .00103 | .00069 | .00046 | .00031 
55 | (00223 | .00409 | .00333 | .00237 | .00168 | .00116 | .00080 | .00055 
60 00493 | .00615 | .00486 | .00355 | .00259 | .00185 | .00132 | .00094 
65 | 00830 | .00865 | .00678 | .00511 | .00382 | .00282 | .00208 | .00152 
70 ‘01225 | .01161 | .00913 | .00710 | .00544 | .00414 | .00313 | .00235 
75 ‘01673 | .01505 | .01197 | .00958 | .00753 | .00587 | .00455 | .00351 
‘80 ~— | .00356 | .02171 | .01900 | .01534 | .01263 | .01015 , .00809 | .00642  .00508 
"35 | 01302 | .02717 | .02348 | .01932 | .01631 | .01338 | .01089 | .00883 | .00714 
"90 ~—s«|-: .02257 | .03310 | .02851 | .02403 | .02068 | .01729 | .01436 | .01188 | .00980 
95 03223 | .03949 | .03412 | .02957 | .02579 | .02196 | .01858 | .01565 | .01316 
1.00 ‘04199 | 04634 | .04035 | .03598 | .03171 | .02745 | .02363 | .02025 | .01733 
1.05 "05186 | .05364 | .04728 | .04325 | .03849 | .03384 | .02959 | .02578 | .02241 
1.10 "06184 | .06140 | .05500 | .05137 | .04618 | .04120 | .03655 | .03232 | .02852 
1.15 ‘07194 | .06963 | .06361 | .06036 | .05482 | .04957 | .04458 | .03997 | .03577 
1.20 (07323 | .07020 | .06445 | .05901 | .05375 | .048S2 | .04425 
1.25 .06956 | .06412 | .05894 | .05407 
1.30 | .07040 | .06531 
ee ET a a 
ing 15 20 25 30 40 50 60 
iis | | 
35 | .00001 | 
40 | .00002 | 
45 =| .00004 | 
50 .00009 | .00001 | 
55 | .00018 | .00002 
60 00033 | .00005 | .00001 | 
65 =| +~-.00059 | .00012 | .00002 | 
.70 00100 | .00024 | .00005 | .00001 | 
75 00161 | .00044 | .00011 | .00003 | 
80 00250 | .00076 | .00023 | .00007 | .00001 
85 00375 | .00127 | .00044 | .00015 | .00002 | 
.90 00347 | .00206 | .00079 | .00030 | .00004 | .00001 
95 00778 | .00323 | .00135 | .00057 | .00010 | .00002 
1.00 ‘01079 | .00489 | .00222 | .00102 | .00022 | .00005 .00001 
1.05 01465 | .00720 | .00355 | .00176 | .0004+ | .00012 .00003 
1.10 ‘01950 | .01033 | .00550 | .00294 | .00085 | .00026 .00008 
1.15 02550 | .O1448 | .00826 | .00474 | .00158 | .00054 700019 
1.20 | .03280 | .01986 | .01208 | .00738 | .00280 | .00108 .00043 
1.25 ‘04155 | .02670 | .01723 | .O1117 | .00476 | .00206 .00092 
1.30 ‘05189 | .03524 | .02402 | .01644 | .00780 | .00376 | .00IS5 
1.35 | .06396 | .04571 | .03276 | .02357 | .01235 | .00656 .00355 
1.40 | .07787 | .05834 | .04379 | .03298 | .01892 | .01098 .00649 
1.45 07333 | .05743 | .04511 | .02810 | .01769 01133 
1.50 07398 | .06038 | .04055 | .02750 .01893 
1.55 | 07920 | .05696 | .04131 |  .03034 
1.60 | | | | .07797 | .06006 .04675 
1.65 | | | | .08465 .06942 
1.70 | | | | | 09949 





2 
Values of k for which P (= < 3) = 0 
8 





n k n k 

4 7811 15 .0468 
5 4775 20 .0259 
6 .3215 25 .0164 
7 .2311 30 .0113 
8 .1740 40 -0063 
9 .1357 50 .0040 

| 10 1088 60 0028 

ll -O891 

12 .0743 
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It will be noted that this function exhibits a more or less singular behavior in 
n — 1 points, the points 


s/s? = at (1 — cos =). 
n 


o~ 5 
p=l,---,n-—1], 


which is particularly well marked for the smaller values of n. 

For odd values of n these singularities are of the order $(n — 4) (i.e. the 
3(n — 3) derivative becomes infinite of order 3). This was shown loc. cit.’, 
p. 390. 

Thus for n = 3 the function has infinities of order $; for n = 5 the direction 
and for n = 7 the curvature have infinities of the same order. 

For even values of n these singularities are also of order 3(n — 4) (i-e., the 
34(n — 4) derivative has an ordinary discontinuity) when uy is odd, but a logarith- 
mic factor must be added to this (i.e. the 3(n — 4) derivative becomes logarith- 
mically infinite) when yp is even. Proofs of these statements will appear later. 

Thus for n = 4 the function has an ordinary discontinuity at ~» = 1, 3 and 
a logarithmic infinity at » = 2; for n = 6 the direction has corresponding 
singularities at 1» = 1, 3, 5 and at uw = 2, 4 respectively. 

For n = 8 (both even and odd) these phenomena would probably be much 
less easily recognized by mere inspection. 








CUMULATIVE FREQUENCY FUNCTIONS! 


By Irvinc W. Burr 


Purdue University 


1. Introduction. The traditional attack upon the problem of determining 
theoretical probabilities and expected frequencies has been through the use of 
the ordinary frequency function. Many such functions have been developed 
for a wide variety of empirical and theoretical situations. The usual procedure 
is to find the “best” function of an appropriate type, and then to integrate 
(either infinitesimal or finite calculus) for the probabilities over the given class 
intervals or other ranges. 

The cumulative frequency function would seem to be theoretically much better 
adapted to this problem. By definition the cumulative frequency function 
gives the expected number of cases less than a given value. Hence expected 
frequencies in any given range are found simply by taking the difference between 
two values of this function. On the other hand, once the ordinary frequency 
function has been determined, these expected frequencies must still be obtained 
by an often times difficult integration. The aim of this paper is to make a 
contribution toward the direct use of the cumulative function so as to utilize 
this theoretical advantage. 

Some properties and theory of the cumulative function will be presented and 
the problem of fitting the function considered. A new cumulative function 
possessing considerable practicability will be discussed and examples given. 


2. Characteristics of the cumulative function F(x). Let F(xo) be the prob- 
ability that « < x. Since probabilities are non-negative, F(x) is non-decreasing 
from F(—«<) = Oto F(«~) = 1. The two ordinary cases will be considered: 
(1). F(x) continuous in (— «, ©) and with F’(x) continuous except for a de- 
numerable set of points in (— «, «), (2). F(x) a step-function with all itsdis- 
continuities at the points nh + d,h >d >0,n = ---, —2, —1,0,1,2,---. 

It is assumed that F(x) has high contact with its asymptotes. Specifically 
for some j (commonly 3 in practice), there is to exist a k > j + 1 such that 
F(x)-a* and [1 — F(z)]z* are ultimately bounded as x tends to — ~ and + 
respectively. These conditions are obviously satisfied when, as is often con- 
venient, a particular expression is used for F(x) over a range, bounded in one or 
both directions, while F(x) is defined =0 below, or =1 above such finite lower 
or upper limits. 

For the continuous case (1), the definition gives at once 


(1) Pia < x < b) = F(b) — F(a). 


1 Presented at the joint meeting of the Institute of Mathematical Statistics, the American 
Mathematical Society and the Econometric Society at Chicago, Ill., September 2, 1941. 
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Furthermore it may be shown that 
() F(a) = [_F(S) 48, F(z) = foo), 
where f(x) is the ordinary probability function. Also 
(3) Paasa<b)=f fle) ae. 


Similarly for the discrete case, 


(4) FQ) = Daf, AFG) = 5f0, 
(5) P(a << b) = Fb +h) — F@) = Dafsid, 


where a, b are among the values nh + d,and A is the usual h-difference. In both 
h 


cases the percentiles are given by the solutions of the equations 
(6) 


Equations (1), (3) and (5) formulate the advantage to the direct use of F(z). 
which was mentioned in section 1. Related to this is the fact that the process 
of finding f(z) from F(z) is at least theoretically much simpler than conversely, 
as (2) and (4) show. The directness of equation (6) is often an advantage also. 

The main problems confronting one in trying to utilize these advantages are 
(a) to find suitable cumulative functions and (b) to find methods of fitting F(x) 
directly. These are next discussed. 






F(x) = n/100. 








3. Some special functions F(x). An obvious method of attack is to use (2) 
or (4) on some f(x). The integration involved is precisely the difficulty the 
writer wishes to avoid. The cumulative function might be sought directly in 
probability theory. A differential equation incorporating some of the properties 
of F(x) given in section 2 is 





d 
(7) = =y(1— y)g(z,y), y= F(x), 
° zx - 
where g(z, y) is to be positive for 0 < y < 1 and z in the range over which the 
solution is to be used. It is to be noted that (7) is very similar to the differential 
equation 


dy _ a = 
- = y(m x)g(x, y), y f(x), 





which generates the Pearson system if g(x, y) = (a + bx + cz’) 






















Ss 
y; 


re 


x) 


2) 
he 
in 
es 


he 
al 





CUMULATIVE FREQUENCY FUNCTIONS 217 


Equation (7) implies the non-decreasing property for F(x), while for many 
choices of g(x, y), dy/dx will be zero at y = Oandy = 1. Wheng(z, y) = g(z), 
(7) becomes 


(8) F(z) = [e470 ® + iy". 


Some functions g(x) whose integrals are such that F(x) increases from 0 to 1 
on the interval —~o < x < @ arec, cx’, [(c — x)z]’, c sec’x and ¢ cosh z, 
where c > 0. Generalizations of their corresponding F(x) are given below in 
(10)-(14) respectively. 

Another method of attack is to simply consider functions which have the 
properties given in section 2. The assumption of high contact provides for the 
existence of certain integrals to be discussed in section 5. Many functions 
having the required properties are to be found in tables of definite integrals, 
particularly Bierens de Haan [1]. 

A list of particular F(z) is given below. In all cases the number of parameters 
would be increased by two by letting x = yz’ + 6, where y and 4 fix the origin 
and scale. These parameters are determined by Z and o. The range of x 
over which the given expression is to be used is written to the right when it is 
not (—«, «). Constants k, r and c are positive real numbers. 





(9) F(z) =2, (0,1), 
(10) F(z) = (* + 1)", 
(11) F(z) = («*+1)", (,-), 
c-—2 Ife —r 
(12) F(z) = ( - ) + 1] , 0), 
—tanz —r Tt 
(14) F(z) = (ke sinhz + a. 
(15) F(z) = 2°(1 + tanh 2)’, 
(16) F(z) = 2 arc tan ¢) : 
2 
ai PO =) nae = +2" 
(18) F(z) =(1-—e*)’, (0, &), 
(19) F(z) = (« -- = sin Dez) (0, 1), 
(20) F(z) =1-(1+2)*, (@,©), 


Most of these functions have unimodal probability functions f(z), and all of 
the functions may be readily handled from the calculational standpoint. To 
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check upon their suitability for practical work, the values of a; and a, for some 
special cases were obtained approximately by evaluating F(x) at a convenient 
regular interval, differencing, and using the results as frequencies of a discrete 





TABLE I 


Calculated a; and a; for special functions F(z) 





| 


Function Parameters a, a, 
(15) r= 1 0 4.01 
(16) r= 1 0 3.24 
(17) k=1,r=2 | — .62 4.50 
(17) k=2,r=1 0 4.11 
(17) k=u2r= 2 — .54 4.22 
(18)? r=1 .63 3.25 
(19)? r=1 0 2.41 


| 
| 
| 
| 
| 





variable. No correction for grouping was made. The values of a3 and a, 
for several of the above functions are given in Table I, where 


[ a’ f(x) dx, +.  f(i) 


t=— 0 


/ 


Kj 


[ @-ud'te@ a, Dae - ws) 


t=— 8 
















a ag F Go = oes 
Co 





It will be seen that a variety of values of a, appear. The values of a3 vary 
considerably in most cases as r varies. These functions show promise of being 
useful after further investigation. The values of a3 and a, for (20) are con- 
venient and adaptable. This function will be discussed in detail in section 6. 





















4. Methods of fitting F(x). The problem of graduation of data by a cumula- 
tive function involves three steps: (a) the selection of the type of function 
(b) the determination of the parameters of the function, and (c) the graduation. 
The first two are often determined by such moment characteristics as a3; and 
a,, as in the Pearson system of frequency functions. The third step involves 
integration or summation if f(x) is used, whereas, once F(x) is fitted, all that 
remains to be done is evaluation of the function and differencing. 

To fit F(x) by moments, it must be possible to determine the parameters of 
F(x) from #, o, a3 and ay. The cumulative moments described in the next 
section, when they can be evaluated, will lead to the values of the Z, o, a3 and a, 
for various values of the parameters. If the relations between the parameters 
and the moments are difficult or impossible to obtain, then tables may be con- 
structed and interpolation used. The usual process would be to use the az 





2 The method of moments of section 5 was used for thesé values. 
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and ay tables to determine the primary parameters such as c, k and r in (9)-(20). 
Then for the given values of c, k, r, one computes the corresponding values of 
= and o from their tables, and these are used to obtain the parameters y and 6 
for 2 = yx’ + 6. This procedure is illustrated in section 6. 

Even when the cumulative moments cannot be evaluated, this method is 
still possible. Graduation by a small interval is used to construct tables of 
¥, o, a3 and a, for varying values of the parameters. Then the table can be 
used as described above. Thus it is seen that in practice any F(x) can be fitted 
by this technique. 

The usefulness of a cumulative or a probability function depends upon how 
wide a range of sets of values of the a; the function covers, and whether such 
values occur in practice. In most of the functions (9)—(20), a3 and a, are con- 
tinuous functions of the parameters. If there is only one parameter then only 
a3 (Or a4) can be fitted in the range of values of a3 which the function possesses, 
but in the case of two parameters both a3; and a, can be fitted. Three or more 
parameters permit a; etc. to be fitted. 


5. Cumulative moment theory for F(z). A moment definition for F(z) isnow 


2 


b 
presented. Since for n > 0, lim / 2"F(z)dzr = =, [ x"F(x) dx cannot be 
bo Ya _ 


used. However, it was assumed in section 2 that for some k > j + 1, 
(1 — F(zx)]z* is ultimately bounded. Hence, lim [1 — F(zx)]z? = 0. Thus 


1 — F(x) can be used as a factor when integrating over any interval (a, ~), 
a being finite. But the factor F(z) must be used for an interval of the type 


(— 2, b). Two integrals are needed, and we define the cumulative moment, 
M ,(a), by 


(22) MU;(a) = / (x — a){l — F(x)] dx — / (x — a)’ F(z) dz, 

which exists under the assumptions of section 2._ The difference of the integrals 
is used because, as will be shown, this leads to simpler results than could be 
obtained by addition. If a = 0 in (22) then calling 17,(0) = M;, 

0 


(23) M; = [ 2’ (1 — F(x)| dx — [ x’ F(x) dz. 


Definitions for the discrete case are similar: 


(24) Mda)=h Dy G— af. — FH) — bh Ds G— FO, 


t=ath t=—00 


0 0 
(25) M; = hd i (1 — FO) — h Don iF, 
t=h i 


i=—o 


where i°”* = i(i — h) +--+ (i —j-— 1h). This function is used because it has 
simpler properties in the finite calculus than has 7’. 
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Various relations between the cumulative moments M ,a) and M;, and be- 
tween these and y;, 4; and a; of (21) are now developed. To express M ;(a) 
7 


in terms of M;’s, use (x — a)’ = >> jC,x**(—a)*. Thus, 
i=0 


M,(a) = [ (2 — a)’{1 — F(x)] dx — [ (x — a)’ F(x) dx 


= . 0 ° a . 
“ [ (fe = ll — Fla — [ (2 = o)' Pie — I (2 — a) dx 
0 a) 0 


= 





ie en Stetina 8 
(26) M;(a) = 2d Ci a)’ Mj; + or: 


One reason for the minus sign of (22) may be noted here, because in the contrary 
case the last term would be [ (x — a){2F(x) — 1)dz. By translating the 


origin in (26) to x = a, renaming the moments, and replacing —a by a, one 
obtains 


qa 


(27) M; = 2 C;a° M,_(a) + —— .€ i 
To bring in ordinary moments, integration-by-parts and (2) are used. 
M,(a) = = {1 — Fa | + Li ea f(a) dz 
(x — a)" ae 


a i [ (x — a)**f(2) dz, 


the first and third quantities vanishing because of the contact assumption. 
A second justification of the minus sign of (22) appears here, since if a positive 
sign were used, the fourth term would have been subtracted and the integrals 
would not combine into (28). Expansion of (x — a)’*’ in powers of x and 
xz — p; yields respectively 


1 j+1 
(29) M,(a) = j+i a iiC(— a)' lieitnts 
1 ’ ; 
(30) M,(a) = i+ rs. 44 ~ iC i(ur _ a)" pips . 
Also setting a = 0, 
1 
(31) M c= mist 


jt+1 

1 7+1 
32 M; = 
(32) 254 


5 oe intCs mn’ Misti 
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It may be shown that the existence of Ma) implies that of the yu; 
i = 1,---,j + 1, and conversely if yu; exists then so do the M,(a) 
+ = 0, ,j — 1. The following formulas are obtained by the opposite inte- 


gration by parts, taking two different forms for | f(x) dx: F(x) and —[1 — F(z)], 
to avoid indeterminate situations. 


, 


Mj f x’ f(x) dx + [ a’ f(x) dx 
—[x’{1 — F(z) }le 
+j [ a’ [1 — F(x)] dx + [x F(a)’. — Jj [ a’ F(x) dz. 


The first and third terms vanish by the contact assumption. Then using 
(c — a + a)*” for 2’, 


j—1 , ‘ 
(33) hj = id jC;a’ M;-1-,(a) +> a, j > 0. 
Also in the same manner 
j-1 a iene 
Mi = FD saCa — wi) Mira) + (@- my’, 5 >1, 


or 


31) my =F D iCd— Mol Myra) + [-Mol@¥, 5 > 1, 


using (29) M,(a) = mm — @. Letting a = 0, 
(35) uj=jMin, j>0 


i-l " > 
(36) Mj=J a 5-C(— Mo)* Mj + (— Mo)’, j> 1. 


An interesting graphical property of F(x) may be seen from (35) j = 1 by 
« 0 
taking »{ = 0. Then M, = 0 and hence I [l — F(s))dr = [ F(x) dz. 


Thus the mean is that ordinate which equates the two areas bounded by (i) 
y = F(z), y = Oand z = yw; and (ii) y = F(x), y = landz = y;. 

It is worth noting that the expressions (34) and (36) have the same coefficients, 
independent of a. This is to be expected because of the invariance of yx; under 
translation. 

If a = »; then (30) simplifies to M (ui) = j + jest: Lastly, expressions 
for a,s in terms of the M;(a)’s are given. 

_ 3M,(a) — 6Mi(a)Mo(a) + 2.Mo(a) 
[2M,(a) — Mi(a)P" 
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The writer has verified that under certain fairly general conditions the dis- 
crete case (38)—(50) approaches the continuous case (26)-(36) as h — 0. 
The following three propositions are merely stated without proof since they 
follow so immediately from (23), (25), (31), (45), (21), (2) and (4). 
PROPOSITION 1: Given a set of functions F(x) and positive constants 


k;i=1,--+,nfor which >> k; = 1, then for F(x) = DO kF x), M, = DOkiiM; 
t=1 t=] t=] 
if all the latter exist. 


Proposition 2. In the above notation, if all the ,u; are equal, then uw; = >> kiiu;, 
t=1 


when the latter exist. 


Proposition 3. If in addition to the above hypotheses, all the :u2 are equal, 
then 


(51) a; = 7 ky sax; - 
t=1 


These propositions are sometimes convenient in forming a linear combination 
of functions F(x), to obtain a function with desired properties. It may be noted 
that Proposition 1 is still algebraically true even with negative k,’s, but these 
might give negative derivatives f(x) for F(x). 

1 
"a 
ulative function will be discussed in detail. The a, can be calculated directly by 
the application of (23), (36) and (21). The resulting a3 and a, values cover a 
broad range, within which those of many empirical and theoretical distributions 
lie. A method of finding such cumulative functions with desired a3; and a, 
will be given. Several graduations are presented for illustration. 

This function appears in Bierens de Haan [1] and has the desired properties. 
The writer has not yet found a probability justification for the function. How- 
ever, since the a; are so close to those of functions which can be so supported, 
it seems that it may eventually prove to be at least some definite approximation 
to a probability situation. 

The complete definition is 


6. An algebraic function, F(x) = 1 This simple algebraic cum- 


1 

fe «i-~ .—. >0 

(52) ” ata = 
= 0 x<Q, 


where c, k > 1 are real numbers. The probability function 


, a ee kexe* 
(53) F(z) = f(x) = G+ 2)” 


is unimodal at x = 4 
ariaeis ck + 1 


l/c 
) ifc > 1, and L-shaped if c = 1. 
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Use of (23) on (52) gives 
° 2g dx ‘ 
(54) M; = | (i + ay’ j<ck-—-—1. 
But from Bierens de Haan [1] 


2 g—l [k—1}¢ 
(55) [ ge dx (c — g)*"* x 


Gite" i-Disswo *** 


where a'« = a(a +c) --- (a+r — le). 
Hence 


(c a j a 1) Me 


ek — 1)!sin2* 


(56) M; 9 ’ J <e-l. 


1 
cs 


However, if j > c — 1 then (55) can still be used through reducing the exponent 
of x by x’ * (1 + 2°) — 2’ * = x’. (56) is only good for integral values of k. 
A more general formula is obtainable by letting (1 + 2°) = 1/s. Then 


1 * j+1)/e—1 k— 1) /e—1 
M; ae Toe (1 on oom 8 (3+) /e ds 
Cc 0 


~ in(it!,,-t4) 
c c c 
1 


J 
rr) 
(57) i a a 


cI'(k) _ 
for j = 0,1, --- up through j < ck — 1, and c, k any real numbers > 1. To 
determine the yu; values the easiest way is to compute the values of the M; 
by (56) or (57), and then to use (36): 


we = 2M, — M, ws = 3M, — 6M\M) + 2M;, 
ie = 4M; = 12M.M, + 12MM; = 3M5 ’ etc. 


Having these, definitions (21) are used for the a;. 

The results for some integral values of k and c are given in Tables II and III. 
These computations were made from (56). Formula (57) shows that for a fixed 
c, M; for k + 1 is obtained by multiplying M; for k by wa = . This re- 
cursion relation is very helpful in the computation, because it enables all of the 
values of the M ,’s for a given c to be found from those for the lowest value of k 
for which M; exists. The values which need to be copied down in the com- 
putation for u;, ¢, a3, a4, by a calculating machine are My, Mi, M2, Ms, 
Mo, M3, Mo, 6Mo, 12Mo, 12M¢, wz, 0, 0, 0, ws, 3, us, a4. Because of 
heavy cancellation, especially in yu; and yu, , it seemed advisable to use eight signi- 
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ficant figures throughout. Eight-place sines were obtained from Gifford [3]. 
The values of the 1/; for k = 11 were also checked by eight-place logarithms 
[4]. These verify the values of the M; for k < 11 because of the recurrence cal- 
culation. 


TABLE II 


anal 
(1 + z*)* 
(In each cell the upper number is a, and the lower number is a) 
| | 


aad 


Mean Pp and Standard Deviation o for F(x) = 1 — 


a 
— |1.20920)1.11072'1 .06896 1.04720 |1.03438 (1.02617 |1.02060 [1.01664 


— | .97787 -58060) .42265-| .33552 | .27953 | .24019 | .21089 | .18815 


\1.00000 | .78540 | .80613| .83304, .85517 | .87266 | .88661 | .89790 .90720 | -91498 
i -= | -61899 | .39533| .30239) .24794 .21116 | .18433 .16375 | .14742 | .13411 


ee | | 








500000) .58905-| 67178) .72891| .769657| .79994 | .82328 | .84178 | .85680 | .86923 


_— .39118 | .29349' .24029 .20461 | 17852 | .15847 | .14253 | .12953 | .11872 


| 33333 | .49087 | .59714) .66817, .71834 | .75550 | .78408 | .80671 | .82507 | .84025 


| 47140 | 30393 | 24784 21077. .18344 | .16234 | .14555-) .13187 | .12053 | .11097 
42951 | .54737| .62641, .68242 | .72402 | .75607 | .78150-| .80215-, .81925- 
| 25596 | .22070, .19269' .17028 | .15220 | .13743 | .12517 | .11487 | .10610 


| | 
| 








| .20000 | .38656 | .51088, .59509)| . .69989 | .73447 | 76196 | .78432 | .80286 
| .24495-) .22488 | .20220| .18010! . | .14505-| .13168 | .12043 | .11087 | .10266 
.16667 | .35435-| .48250) 57029 68045-| .71698 | .74609 | .76980 | .78948 
.19720 | .20274 | .18851| .17064) . .13962 | .12731 | .11682 | .10782 | .10004 








| -14286 | .32904 45952) .54992) . .66425-| .70235 | .73276 | .75758 | .77820 
8 | .16496 | .18599 | 17783) .16316) .14846 | .13528 | .12383 | .11394 | .10539 | .09796 
| .12500 | .30847 | .44038) .53274) .59982 | .65041 | .68981 | .72131 | .74706 | .76848 

.14174 | .17275 | .16918, .15704' .14389 | .13171 | .12095-| .11156 | .10338 | .09623 
| | | | 


' nes i Ce ee ee eee eT ee ee 
| .11111 | .29134 | .42407| .51794 .58649 | .63836 | .67886 | .71130 | .73783 | .75994 
| -12423 | .16197 16197; .15190: .14002 | .12869 | .11851 | .10954 | .10168 | .09477 























.10000 | .27677 | .40993 .50499| .57476 | .62772 | .66916 | .70240 | .72964 | .75234 
11055 | .15297 | .15584 .14749 .13669 | .12608 | .11640 | .10780 | 10021 .09351 








It will be seen from Table II that in most cases the values of a3 and ay lie 
within useful ranges. The graph shows the general relationship between a; , 
k and c. The curves are the traces of planes k = 1, 2, --- upon the surface 
a; = G(c, k). Other traces would contain all pairs (c, k) giving a fixed a;. 
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TABLE III 
seine 
(1 + x°)* 


(In each cell the upper number is a; and the lower number is a4) 


Skewness a3 and Kurtosis a4 for F(z) = 1 — 





1.820) 1.458 


190 | 












. 207 


4 —_ 7.356 | 4.036 | 3.363 | 3.189 | 3.169) 3.205 | 3.263 | 3.327 | 3.393 


.648| 1.218 | .559 | .238| .040 | —.097| —.199 | —.277 | —.340 | —.391 
80 | 5.832 | 3.604 | 3.154 | 3.070] 3.098} 3.165-| 3.243 | 3.324] 3.401 
SC ES a eS 
—.013 | —.147| —.246 | —.323 | —.384 | —.435 
6 | 38.67 | 5.118 | 3.380 | 3.045 | 3.010 | 3.065) 3.150-| 3.241 | 3.330 | 3.416 


| 
































hentia gees -avnafhstonacrailsodesnanliniones 
3.381) 1.014 .433 .136 | —.051 | —.181) —.279 | —.355-| —.415 | —.465 
7 | 27.86 | 4.707 | 3.245-| 2.979 2.975 | 3.048) 3.144 3.244 | 3.339 3.430 


| 
| 





| 3.118 .958 | .396 | .106 | —.078 | —.207; —.303 | —.378 | —.438 | —.488 


8 | 22.73 4.443 oe 2.953 | 3.039) 3.143 | 3.248 | 3.349 | 3.442 
























































ees - eae cer recs 
| 2.940) .916 368 .083 | —.098 | —.226) —.322 | —.396 | —.456 | —.505- 

9 | 19.76 | 4.258 | 3.091 | 2.906 | 2.938 | 3.033) 3.143 | 3.252 | 3.357 | 3.453 

| 2.811) .884 347 -065 | —.115-| —.242) —.336 | —.410 | —.470 | —.519 

10 | 17.83 | 4.122 | 3.043 | 2.883 | 2.928 | 3.030) 3.144! 3.257 | 3.364 | 3.462 




















| 2.714] .858 om 


329 0580 | — .254| —.348 | —.422 | —.481 | —.530 
11 | 16.48 | 4.018 | 3.006 | 2.866 | 


128 
.920 | 3.027) 3.146 | 3.261 | 3.371 | 3.470 


to 














, . . e 
The surfaces for yn; , o and a, are more irregular. The problem of determining 
a cumulative function with a; = a and a, = b is equivalent to the problem of 
determining a point of intersection of the curves 









a3 = G(k, c), a3 =a 


a, = H(k, c), a, = b. 


(58) 


Direct algebraic solution of this system appears verv difficult, and other tech- 
niques must be resorted to. 
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One method is to use only integral values of k, and then for each k interpolate 
for the value of c giving the desired a3. For such pairs of c and k, find a, by 
interpolation. Then choosing the pairs having a, just above and just below the 
desired one, the proper linear combination (51) is taken. This gives a combina- 
tion function which has both a; and ay at the desired values. This combination 


Figure I 


will be an approximation to the single function with non-integral k, having the 
given a; and o,. This method of linear combinations might be extended to fit 
as by using three integral values of k. 

The interpolations may be done graphically by use of Figure 1 and others like 
it. Or one may use Stirling’s formula [5]. The interpolation for c from ag 
is backwards, while that for a, from c is direct. Sometimes it is more accurate 
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to use Newton’s formula [5, p. 36] when the values in one direction increase 
rapidly. 

Use of a single function F(x) for a graduation is easily accomplished. First, 
obtain the c and k to be used so that az is correct and a, is as close to the given 
value as possible. Then determine y; and o from Table II by interpolation. 
Change the scale and origin of the original values of the variable X to those 
x’s corresponding to F(x) = 1 — 1/(1 + 2‘)*, through 
(59) aah ee 


ine See ae 


o i. 








where M and S are the mean and standard deviation of the given distribution. 
Now compute the values of 1/(1 + 2°)“ for the various values of x. The differ- 
ences of these results are equal to the differences of F(x), which by (1) are the 
probabilities for the given ranges of X. Multiplication by the total frequency 
will yield the theoretical frequencies, if desired. 

If the graduation is to be done by a combination of two functions, the work 
is carried out for each as described above, and then the frequencies are combined 
by the same linear combination as that by which the component a,4’s must be 
combined to give the desired a,. This may readily be seen by considering the 
separate cumulative functions in terms of the standard variable ¢, whence the 
means and o’s are 0 and 1 and (51) is applicable. Then the differences of 
G(t) = kiG,(t) + keG.(t) are sought. But these can be found by taking the same 
linear combination of the separate differences of the functions G,(t) and G,(t). 
However, these values are merely computed from their respective sets of x values. 

For illustration, three graduations are given. The first is a highly normal 
distribution of heights from Rietz [5, p. 98ff.]. For this distribution, M = .02085, 
S = 2.5723, a3 = —.0124, ag = 3.149. The graduation was done by taking the 


function F(x) = 1 — see seee which has the nearly normal characteristics 
(1 + 2*)* 
a3 = —.019, ag = 3.169. The object was to take a simple cumulative function 


with integral k and ¢ to show how a satisfactory job can be done on a normal 
distribution. For this function n; = .75550 and o = .16234. Then 


x = .063110X + .75418, 


into which are substituted the X class-limits —11.5, —10.5, ete. 
(i ee , are calculated and differenced to give the 
theoretical frequencies for the 8585 cases. The results are given in Table IV. 

The fit obtained by use of F(x) is good. One comparison test is that of 
x’. The eight classes —11, —10, —9, 9, 10, 11, 12, 13 were grouped together. 
The results were 


From these, 





corresponding values of 





x; = 21.210, xy = 23.479, 
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as compared to 
P(x’ > 22.31) = .10, 
P(x? > 19.31) = .20, 


for 15 degrees of freedom (18 classes minus 3 for linear restrictions). One 
reason for the somewhat lower x’ for F(x) may be that its a3 and ay are closer to 


TABLE IV 





Observed frequency [5] Cute by | eee aed of 





.00 | 168 
43 | .67 
2.84 
10.30 
32.11 


86.03 
198.17 
392.43 
668.11 
977 .92 


— ry 


1230.63 
1331.41 
1238.41 
990.33 
680.86 


bl ss we Ay Vw 





402.44 
204.51 
89.35 
33.56 


$585.00 8584.99 


those of the observed distribution than are those of the normal function. This 
gives a better fit in the tails of the distribution. Nevertheless, this example does 
illustrate how one of the simplest of the cumulative functions with “normal” 
characteristics can be used without specifically fitting a; and a,;. It may also be 
mentioned that F(x) for c = 5, k = 6 has a3 and ay even closer to the normal 


3 Total of stump frequency. 
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Observed 
frequency i 


3 .00 

9 | 86 
46 39.58 
167 180.78 


1186 1116.06 
1462 1383 .06 
1498 1492.04 


1419.70 





1205.81 
913 | 926.59 
642 | 654.00 
435 430.66 


268.70 





161.10 


133 93.99 
47 | 53.88 
29 30.62 


17.37 


9.86 
5.64 
3.26 
| 1.89 
1.12 


ND © oro 


| 66 
| Al 
| 24 
| 16 
27 
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Type III [6] 




















TABLE VI 


Type A [6] 


i) 
Roo 
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2 
os ff 
~ 
o 


| F(z) = .3063/. 


+ .6937fs Type IIE [5] 








00 
25.07 
175.27 
445.79 


791.72 
1134.52 
1384.99 
1477.86 
1399.70 


1190.80 


921.47 
656. 82 
436.78 
274.63 


165.27 


96.23 
54.77 
30.70 
17.07 


9.46 
5.26 
2.93 
1.66 
.93 
.53 
31 
18 
aa 
174 


| 
| 














Edgeworth [6] 







| .00 

wT i 2 

| 29.52 7 

| 176.96 142 
442.13 











784.71 


1128.86 1186 
1384.40 1441 
1482.20 1502 


1405.83 







1195.40 


923.01 891 
655.96 641 
434.90 434 


272.81 


163.99 | 


95.55 | 102 
54.50 | 59 
30.68 33 


17.16 













9.58 9 
5.38 5 
3.03 2 
1.73 | 1 

.99 1 





57 | 
34 | 
20 | 
13 | 







10700 .99 





10701 






F(x) 


59 56 
53 52 
32 34 
15 16 








‘Stump frequency. 
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values, but it does not give quite as good a fit because it tends to decrease too 
rapidly on the left. 

The second example is also from Rietz [5, p. 108ff.]. For this distribution, 
M = .68835, S = 2.9480, a3 = .583 and a, = 3.698. Two functions were used 
with k = 4andk = 5. By interpolation 


My o a3 as 
k=5 c = 2.944 .54200 .22247 583 3.655 
k=4 c = 3.228 .61577 .23823 583 3.795 


Because of the rather rapid increases for smaller values of c, Newton’s formula 
[5, p. 36] yields better approximations than Stirling’s [5, p. 38 (12)]. The gradua- 
tion for each function is carried out as above, and since 


.3063 -3.795 + .6937-3.655 = 3.698, 
the linear form 


.3063f; + .6937f; = f", 
is used. 

Table V gives the component and combined frequencies, and also the fre- 
quencies froma Type III.  ° for both are very high even though the fit appears 
reasonably good ona graph. This result is due to classes 6 and 8 which tend to 
cause a high x’ for any distribution function of a small number of parameters. 
The example, however, does show that F(x) can be used to graduate a skewed 
distribution. 

It is to be further noted that the component functions were used only to 
obtain an approximation to a single function with + < k < 5, for which a; and a, 
are simultaneously correct. When tables more complete than Tables II and 
III are available, such a single function can be found. 

The third example of graduations is from Elderton [6]. The measures were 
treated as a discrete variable in computing a3 and a,. A single function 
c = 3.102, k = 11 was used. This function had a; at the observed value of 
.2936, while a; was 2.973 as compared to the observed 2.986. The results along 
with those by classical methods are shown in Table VI. The above x’ were 
obtained by grouping the first and the last three class frequencies. The values 
are approximate because of rounding. However, they do show that F(x) does 
a comparable graduation. 

Besides aiding in the problem of graduation, this cumulative function should 
prove of value in the approximation of known or population distributions, as 
for example, (p + q)". However much more work needs to be done before this 
can be more than a conjecture. 


7. Conclusion. This paper has stressed the advantages obtained by the direct 
use of the cumulative function. A number of useful functions have been 
considered. A general method for fitting any cumulative function by the 
construction of a table has been suggested. A particular method depending 
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upon the use of certain new cumulative moments has been given. Making use 


of this theory a certain simple algebraic function has been discussed in detail, 
and its use in graduations explained. 


The writer wishes to convey his sincere thanks to Professor Harry C. Carver, 
whose counsel was most helpful. 
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NOTES 


This section 1s devoted to brief research and expository articles, notes on methodology 
and other short items. 


— gg 


AN APPROXIMATE NORMALIZATION OF THE ANALYSIS OF 
VARIANCE DISTRIBUTION 


By Epwarp Pavutson’ 


Columbia University 


The statistic F = s{/s} , where sj and s3 are two independent estimates of the 
same variance, has played an essential part in modern statistical theory. All 
tests of significance involving the testing of a linear hypothesis, which includes 
the analysis of variance and covariance and multiple regression problems, can 
be reduced to finding the probability integral of the F distribution. This 
distribution (and the equivalent distribution of z = 3} log F) has so far been 
directly tabulated only for the 20, 5, 1, and 0.1 percent levels of significance [1]. 
To find the critical value of F for some other probability level would require 
the use of Pearson’s extensive triple-entry tables [2], which is not very con- 
venient to use for this purpose, and in addition is inadequate for some ranges 
of the parameters. 

It therefore appears that it might be of some practical value to have an 
approximate method of determining the critical values of F for other probability 
levels. A solution will be given based on a modified statistic U, a function of F, 
so selected as to tend to have a nearly normal distribution with zero mean and 
unit variance. This normalized statistic will have the additional advantage 
that further tests are possible with normalized variates, as pointed out by 
Hotelling and Frankel [3]. 

F can be written in the form 


iat xi/m 
x3/M" 


where xj and x3 have the chi-square distribution with n,; and n2 degrees of freedom 


\4 
respectively. It is known from the work of Wilson and Hilferty [4] that (x) 


is nearly normally distributed with mean 1 — 2/9n and varianée 2/9n. An 
obvious approach to the problem of securing an approximation to the F distribu- 
tion is to regard F* as the ratio of two normally distributed variates. In general 
the distribution of the ratio v = y/x where y and z are normally and inde- 
pendently distributed with means m, and m, and standard deviations o, and o, 


1 Work done under a grant-in-aid from the Carnegie Corporation of New York. 
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is not expressible in simple form. However Fieller [5] has shown that a function 
ms will be nearly normally distributed with zero 
Vy? o: + a, 

mean and unit variance, provided the probability of x being negative is small. 
In the given problem it follows that we can regard 


(-2)"- 0-8) 
(1) ty = 9n. : 9n, 
re. 
VES: + 9n, 


as nearly normally distributed (with zero mean and unit variance) provided 
n2 > 3, for with n. = 3 the probability of the denominator of F* being negative 
is only .0003. If it is desired to use the lower tail of the F distribution, then the 
statistic U should only be used if m; is also > 3. Ordinarily, in most applica- 
tions only the upper tail of the F distribution is used, and nz , which corresponds 
to the number of degrees of freedom in the estimate of the error variance, will 
be much greater than 3. 

The following tables show the degree of accuracy of the approximation. The 
exact value of F corresponding to various levels of significance are compared 


R of v, namely R = 





m=1, m=10 
t=/F 


Approximation | Exact Value 


| 
i 
| — 


n, = 4, ne=8 


F 


Approximation Exact Value 





161 

.407 
1.92 
3.84 


.058 
.12 | 
| 


14.39 


1.37 
2.21 
3.16 
4.63 
6.40 


.068 

- 166 

. 406 
1.92 
3.84 
7.01 


1.37 
2.23 
3.17 
4.59 
6.22 


Approximation Exact Value 
.123 
.248 
497 
12 
.00 
.85 
.58 
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with the approximate values, which are found by solving (1) for F by considering 


it as a quadratic equation in F*. In these tables P = i ¢(F) dF, where ¢(F) 
F 


is the probability distribution of F. The case n, = 1 is of special interest, 
since here F = ¢’, where ¢ has Student’s distribution, and is shown separately. 
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NOTE ON THE DISTRIBUTION OF ROOTS OF A POLYNOMIAL WITH 
RANDOM COMPLEX COEFFICIENTS 


By M. A. GrirsHick 
United States Department of Agriculture 


In order to obtain the distribution of roots of a polynomial with random 
complex coefficients, it was found convenient to employ a rather well known 
theorem on complex Jacobians. Since proofs of this theorem are not very 
plentiful in the literature, a brief and simple proof of it is presented in this note: 

THEOREM: Let n analytic functions be defined by 


(1) Wp = Up + Wp = fol%, 22, °** 5 Zn); (p = 1,2,---,n), 


where 2p = tp + yp, i = W—1. Let j denote the Jacobian of the transfor- 
mation of the n complex variables defined by (1). That is 


_ Ow, | 


(2) 


Let furthermore J denote the Jacobian of the transformation of the 2n real variables 
defined by the equations up = Up(%1, T2,°**y Lnj Yry Y2> °°» Yn) aNd vy = 
Up(%i , Ze, eee >In, Yr, Y2,°°° » Yn); (p = 1,2,--- ,2).. That is 

Us Uy 


(3 J=) , 
'V. Vy! 
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Then J equals the square of the modulus of j. 


ow ; 
Proor: Since by hypothesis w, is analytic we can set —? = —? —i —?. 
y yp Pp y d2Zq dYq OY ¢ 
Hence j takes on the form: 













(4) j=|V,—iU,|. 
. , . . OUp _ Wp Wp _ _ IUy hat j 
Again, since w, is analytic, we have in * ce.’ & a That is 
U, = V, and V, = —U,. Hence J in (3) has the value 
(5) ra "|. 
| = U, Vy | 


.Now J can also be written in the form 


| Vy 
(6) J=|. 
\2U, Vy | 





This follows from the fact that if we multiply each of the last n rows of the 
expression for J in (6) by 7 and factor out 7 from the last n columns, we get the 
expression for J given in (5). 

Now in (6) subtract the (n + p)th row from the pth row for each p = 
1, 2,---,n. This yields: 





|V,—iU, iU, — Vy! 
| aU, Vy 


(7) J 










Next add in (7) the pth column to the (n + p)th column for each p = 1,2, --- ,n 
This yields: 





'V—iU 0 


(8) J =| }=|V—-—iv||V+iU I. 
iU V+iU| 
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But (8) is precisely the square of the modulus of |V — iU |. This in con- 
junction with (4) proves the theorem. 
Consider the equation 


(9) 2” — az” +++. + (—1)"a, = 0, 


where the a; are complex numbers. We may wish to consider the real and 
imaginary parts of a; as random variables having a given joint distribution 
function, and require to find the probability that one or more roots of (9) will 
lie in a specified region of the complex plane. In order to answer this question, 
it is necessary to find the joint distribution of the real and imaginary parts of 
the roots of (9). 

As an example let us assume that the real and imaginary parts of a, are nor- 
mally and independently distributed with zero mean and variance o°. That is, 
we assume that the distribution density of these quantities is given by 


1 2n 1 n _ 
(10) (se) exp |- Do? a 254, | 


where 4d, is the conjugate of a,. Let 2, 2, --- ,2, be the roots of (9). The 
relationship between the roots and coefficients of (9) are given by 


nn 
(11) = 223, 2 = De zjeey sss, On = are An 
j=1 1<k 


Thus the a,’s are analytic functions of the z’s. 

In order to find the joint distribution of the real and imaginary parts of the 
z’s, it is necessary to find the real Jacobian J of the transformation defined 
by (11). Now the complex Jacobian j of the transformation (11) is defined as 


jaa 


(12) 


A simple calculation will show that the value of j in (12) is given by 
(13) j= dL L &— a). 
p=1 g=ptl 


Hence, applying the theorem proved above, we get 
(14) J=(|jP= 2 » | zp — zl’, 


p=1 g=ptl 


where the symbol || stands for the modulus. 
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From (10) and (14) we conclude that the joint distribution density of the real 
and imaginary parts of the roots of (9) is given by 


1 2n 1 n n 7 
(Far) op —aa{ee eat 
+ at e63s} | > |e— 2%’. 


p=1 g=ptl 
















(15) 





A NOTE ON THE PROBABILITY OF ARBITRARY EVENTS 


By Hixpa Gerrincer' 
Bryn Mawr College 





In a recently published paper [1] on arbitrary events the author studies the 
probability of the occurrence of at least m among n events. Denoting by 
Pm(¥1, ¥2, *** Yr) the probability that at least m among the r events, E,, , 
--+ H,, occur, and by pja,,a:,...2,) the probability of the non occurrence of the 
events numbered a; , a2, --- a, and of the occurrence of the n — r others, he 
proves 





— PrlOrst » —? Gn) + Daly, Gr+1,°°° On) _ Zo ~ may, Y2, @r41,°°° Qn) 
ve v1 


v2 


(1) 





















+--+ + (-1) D pill, +++ 0) = ptay,---a,1- 


(Theorem VI, page 336). From (I) he deduces that a necessary and sufficient 
condition for the existence of a system of events E,, --- E, associated with 
given values ¢,; (a1 , --- ax) is that the expressions on the left side of (I) computed 
from these ¢’s are 2 0 for all possible combinations of the a’s (Theorem VII). 
He also points out that it was not possible to find similar (necessary and suf- 
ficient) conditions for m ~ 1. I wish to show in this note the relation between 
these theorems and some well known basic facts of the theory of arbitrarily 
linked events and to add some remarks. 


1. Given n chance variables x; (i = 1, --- n) denote by x; = 1 the ‘‘oc- 
currence of E,’’, by xz; = 0 its non occurrence and by v(x, 22, --- 2m) the 
probability of “the result (1, 22, --+ 2%n)’’ ie., the probability that the first 
variable equals x, the second 2, --- the last z,; e.g. v(1, 1, 1, 0, --- 0) = 


Vi4s...n} 1S the probability that only the three first events occur. Hence the 
v’s are 2” probabilities, arbitrary except for the condition to have the sum 1. 

Instead of these v’s we often introduce another set of 2” — 1 probabilities, 
namely p; the probability of the occurrence of E; (¢ = 1, --- n); p;; that of 
the joint occurrence of E; and E; (i,j = 1, --- n); --- pi... the probability 
that all the events occur. 





' Research under a grant-in-aid of the American Philosophical Society. 
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It may be noted that instead of the p;, pij, --+ piz...n we could quite as 
well use a system of q;, Gij, °°* Qiz.--n where qi is the probability of the non- 
occurrence of E; (or of the occurrence of E; = E — E;, qi; that of the joint 
non-occurrence of E,; and E; (of the occurrence of E;E;) and qy..., the probability 
of ExE; ‘se E’. 

The use of the p’s (or q’s) instead of the ‘elementary probabilities” v is 
justified by the fact that the p’s are (2" — 1) independent linear combinations 
of the v’s and that therefore the v’s and the p’s (or the v’s and the q’s) determine 
each other uniquely. There exist in fact the following well known relations, 
(1) and (2). The first set (1) gives just the definition of the 2” — 1 probabilities . 
p, in terms of the v’s, andthe second set expresses the v’s by the p’s as the result 
of the solution of the 2” — 1 independent linear equations (1). Thus we have, 
beginning with py... 


Pi2---n = v(1, 1, <e 1), 
Su... 2d (i, 1, ), 


°° » Do v(a, tay *** Seay), 
and solving successively: 


v(1,1,--- 1) = pp... 
v(1, 1, 


= 2d Prin +> 2X 2X Prin — 


= * ->d Priva-+-Yn—17 - Piz.- 


Yn-1 


The successive solution of the system (1) with respect to the “unknowns” v 
is possible because each new equation in (1) contains exactly one new unknown v; 
e.g. in the equation defining py the only “unknown” is v(1, 1, 0, 0, --- 0) all 
the v’s with more than two “1”’s having already been computed from the fore- 
going equations. 

If we choose to use the system of the q’s we have in the same way: 


7. 
: ** » Do v(t, te, +++ Inn, 0); 


Zam) 
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and the inverse system 


v(0, 0, «++ 0) 


= du. .-n 


(2') v(0, 0, Pe 0, 1) = Gie..-n—1 — QN2.--n 


eee eee eee eee eee eee eee eeeeeseeees 


v(1, 1, a 1, 0) == Eten + 7 he Qyiven — °°° + Qiz...n- 
v1 vs 6S 
Coming back to Chung’s theorem we see that the probability pi(ai, --+ ar) 
that at least one event among E,, , --- Ea, occurs is evidently: 


(3) Pia , ah a,) =l- Jaj,...ap o 


If we introduce this value in (I) all the “‘1”’s introduced by (3) cancel and we 
get our system (2'). (Of course we could in the same way deduce from (2) a 
system of equations for gi(a:, +--+ a,) = 1 — Pa,,...a, Where gi(ai, -** a, is 
the probability that at most r — 1 events among the r given ones occur.) 

As the v-values on the left side of (2’) are 2" — 1 independent probabilities 
only subject to the restriction that they have the sum < 1, we sce that the 
expressions on the right side of (2’) must have the same properties; and these 
properties are also sufficient for a system of q’s, (or for a system of pi(ai, --- 
a,); indeed if they are fulfilled, these 2” — 1 expressions define by means of 
(2’) a system of elementary probabilities v(z, --- x,). Hence the theorems 
VI and VII quoted at the beginning of this note are rather close consequences 
of the basic relations of the theory of arbitrary events. 

2. Remark 1. We may add one more equation to equations (1), namely 


1= Z. wel > v(x, *** Zn), 
Z Zn 


thus introducing v(0, 0, --- 0). Then in system (2) the corresponding new 
equation will be: 


° 


»0,0,---0) =1- hs tn + ) > a — tts  DPr...n 
v1 va 6 OTS 


(and analogously for the q’s). In this way we get two systems (1) and (2) 
each consisting of 2” equations and in (2) the sum of the expressions on the right 
side is now identically equal to one. Hence necessary and sufficient conditions 
will now be that all these 2” expressions must be non-negative. 

Remark 2. It is convenient to interpret or prove results of the kind con- 
sidered here in terms of elementary measure theory: p; is the measure of a set 
E, ; po of Ez; py that of the intersection £,E2 etc., and analogously for the 
v’s: e.g. v(1, 1, 0, --- 0) = m(E,E2E; --- E,). Consider now the equations 
(2). The first is an identity. In the second py....-1 measures the product of 
E\E, +++ Ex, whereas piy...n-1 — Pi...n is the measure of that part of this 
product which does not belong to EF, , and it therefore equals m(E,E, --- 
E,,,E',) = v(1,1---1,0). Inthe last equation (2) >> --- >> Pv, .-:ta~gn 8 the 

v1 


Yn-2 
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measure of that part of Z, which belongs to at least (n — 2) other sets (besides 
E,,); whereas this same value minus pi..., is the measure of the part of E, 
which belongs exactly to (n — 2) other sets; now subtracting this expression 


from >, --: Py;..-%m-3 We get the measure of the part of EZ, which belongs 
Yn 


v1 =% 

exactly to (n — 3) other sets and finally p, — --- + py..., is the measure of 
that part of E, which belongs to no other set besides, i.e. m(E;E; --- E,+E,) = 
v(0, 0, --- 0,1). This kind of proof does not require the solution of (1). 

Remark 3. According to (1) the p;, pi;, +++ pie.... are the ordinary mo- 
ments of order 1, 2 --- n of v(x, 22, +++ tn). There are of course many 
more than 2” — 1 moments of this n-variate distribution but only 2” — 1 of 
them are different from each other because 1” = 1. 

3. Denote by p,(x), (x = 0, 1, --- n) the probability of getting exactly x suc- 
cesses in n trials. (See e.g. [2], [3].) For the simplest case of arbitrary events, 


the Bernoulli problem, p,(x) = (") p (1 — p)”*. Then the probabiliy of 


at least x successes (of a number of successes = 2) is 
(4) Vi(t) = pa(x) + pa(x + 1) + +--+ pa(n), 


or p:(1, 2, --- nm) in Chung’s notation. The p,(x) are by their definition 
(n + 1) arbitrary positive numbers with sum equal to one. These are the only 
necessary and sufficient restrictions for p,(x). V,(x) the “cumulative” distri- 
bution of p,(x) which is defined for x between (— © and +) is a monotone 
non-increasing step function with its (n + 1) steps at x = 0, 1, 2, --- n equal 
to the p,(z). 

Consider next p.(a1, a2, --* a) where r < n; these are cumulative dis- 


tributions each corresponding to one of the (") probabilities p,(z) where p,(x) 

is the probability of exactly x successes in a group of r trials.” For each group 

(a1, °*+ @,) the corresponding p,(x), (x = 0,1, --- r) are positive and with 
O--er 

sum equal to one. Hence if we always omit p,(0) because of >> p,(z) = 1, 


all the different p,(x), po(x), --- pna(x) together define 
1l-n + n(n — 1) + (3) (n—2) + --- +n = n2”™" 


values. As n2"' > 2” for n > 2 we realise that between these n2” prob- 
abilities there must exist a set of n2”-* — (2" — 1) identical relations; and the 
same is true for the corresponding cumulative distributions V,(x) or p.(a, 
- a,). Thus it seems reasonable that it may be hard to use these p.(a, 
- ar) in the characterization of a problem of arbitrarily linked events if 
x > 1. On the other hand we have seen in 1 that for x = 1 they reduce to the 


2? One may write here p,(z) instead of pcaj,a2,°-+,a,)(Z). 
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2” — 1 probabilities q:, qi;, --* Qi2.... Which of course define the system of 
events unequivocally. 


4. Introduce in the usual way the sums of the p; , pi; , ete. 


(5) Si = DP, S: = De Piss --- §, = Pi2.--ny and Sp = 1. 


Now add in system (1) first the n equations which define p,, po, --* Dn, 
then the “ equations for the p;;, etc. Observing that p,(x) is the sum of 


all these elementary probabilities v(x: 22 --- Z,) with exactly xz “‘1”’s and (n — 2) 
“0’’s we get as the result of these n additions the well known formulae: 


(6) S, = > o pit), (y =0,1,--- n). 


z=7 
n 


Here y = 0 gives Sy = 1 = >> p,(x). We may solve successively these (n + 1) 


0 
linear equations with respect to p,(n), p.(n — 1), --- p,(O), each linear equa- 
tion containing only one new unknown, and ‘ind: 


(7) pr(x) = >, (—1)"* (7) S., (x = 0,1, --- n). 
y=2 

(These formulae could also have been derived from (2) by collecting groups of 
equations such that all the corresponding v(x, --:- 2,) contain the same 
number of ‘‘1’’s.) (In the measure interpretation p,(x) is the measure of that 
part which belongs exactly to x of the original sets and S, measures the set which 
belongs to at least 7 of these sets.) We also find by “‘cumulating”’ equations (9) 


(8) Sy ™ ie (? . ') V(x), (y = i, 2, wk n), 
zy \Y — 1 

and the inverse system 

(9) V,(z) = Zz (yr (? _ + S,, (x = l, 2, re n). 
y=z 


(6) and (8) are of the same type as (1), and (7) and (9) of the same as (2). We 
also may deduce analogous formulae by interchanging the roles of 0 aud 1 and 
introducing a system of 7, T2, --- T, which depends on the q’s in the same 
way as the S,, S2, --- S, defined-in (5) depend on the p’s. 

We have seen that the p,(x) are (n + 1) arbitrary nor-negative numbers 
subject to only the condition of having their sum equal to one. But the S, 
(y = 0, --- ”) are not arbitrary as we see from (7). The (nm + 1) expressions on 
the right side of (7) must each be non-negative if they are to define the probabilities 
p,(x) (their sum is identically equal to one). Then and only then they define 
a system of arbitrarily liked events E,, --- E,. 

The p,(x), (cx = 0, 1, --- n) are of course not equivalent to the complete 
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system of 2” — 1 values v(x, , 22, --- 2,) and the same remark holds for the 
So, -:: S, and the system of p;, pi;, --+ Pw...n- But often we are par- 
ticularly interested in problems dealing only with the p,(z) (and S,). (For 
instance the author has studied [2] the asymptotic behavior of p,(x) as n tends 
in different ways towards infinity.) The simplest way to indicate a particular 
p-system corresponding to given S, is of course to assume all the p,; equal to 
each other; all the p;; equal to each other etc. and to put therefore: 


1 
aR Oh om, 


Pr = eco SS Pa-i,2 = [1 /(3) | eee Pi2...n = Sn. 


In the corresponding v-system all these v’s which show the same number of 
“‘1”’s equal each other. 

We see from (6) that the S, (multiplied by y!) are the factorial moments of 
order 0,1, --- ” of the distribution p,(x). Therefore by (7) we get the p,(z) 
in terms of their factorial moments up to order n. We may therefore also say: 
Necessary and sufficient conditions that a system of numbers No = 1, Ni, °°: Na 
be the factorial moment of an arithmetical distribution with at most (n + 1) steps 
atx = 0,1, --- n are the inequalities: 


(10) . Sees enone 


z!(y — 2)! = 

Note that here there is no more allusion to a set of arbitrary events; (10) are 
the necessary and sufficient conditions for a set of (n + 1) numbers to be the 
(n + 1) (factorial) moments of an arbitrary arithmetic distribution with its 
abscissae given. The linear inequalities (10) differ very much from the basjc 
inequalities in the classical problem of moments; because in our problem the 
abscissae of the steps are given in advance. 


5. In some problems (e.g., some questions connected with the law of large 
numbers, with correlation theory, with analysis of variance) we are only con- 
cerned with the first and second moment of a distribution. Thus we are lead 
to the following question: Given r + 1 numbers No, Ni, --: N,, (r S n) 
indicate a set of necessary and sufficient conditions such that these numbers 
are the moments of an arithmetic distribution with at most (n + 1) steps, at 
0, 1, 2, --- n.*> Some sort of an answer which may work well in particular 
cases, can immediately be deduced from (10). “r + 1 numbers No, Ni, 
--+ N, will be the factorial moments of an arithmetic distribution with, at 
most, (n + 1) steps at 0, 1, 2, --- n if and only if it is possible to indicate s 


3.This problem and the method of its solution has much in common with a problem 
studied in R. von Mises’ paper [4]. 
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numbers N,41, --: Nii, (0 S$ s S n — 1), such that for the r + s + 1 num- 
bers No, Ni, --- Nr, +++ Nrys the r + s + 1 inequalities 


r+s (— 1" 


4 a ar ** (x = 0,1, ---r +s) 


be satisfied.” 
The proof of this statement is self evident but the statement itself cannot 
be considered satisfactory. We get a general solution in the following way. 


Let fi(t), --- f-(t) be r functions of the chance variable ¢, v(¢) an arithmetic 
probability with n given attributes t;, f, --- t, and 


(11) El f,(t)] = > fultso(ty) = Ye ay0, = Bos (p = 1, 2, me r), 


the expectations of f,(¢) with respect to v(t). We wish to indicate necessary 
and sufficient conditions for the r numbers S,. For f,(t) = ¢’ we have the prob- 
lem stated above where the first r moments are given. 

Call (S) the r-dimensional curve x, = f,(t) and P;, P2, --- P, the points 
on (S) with coordinates f,(t,) = a,,, (op = 1,---r;y = 1,--- 7), S the given 
point with coordinates S,. In this case, the point S must be contained in the 
smallest convex body (B) determined by the n points P,,--- P,. This condition 
as necessary and sufficient. Because, if we interpret the v, which are 2 0 as 
masses of the points P, , with sum equal to one, then S is the center of gravity 
of these masses and it is well known that the above mentioned condition for S 
has to be fulfilled. But this condition is also sufficient, because if S is contained 
in (B) there exists always a simplex of at most r dimensions, consisting of at 
most (r + 1) of the given points such that S is the center of gravity of appro- 
priate masses in these points. 

If we want to indicate explicitly the inequalities for the S, we must know 
the boundary of (B). This is determined by its planes of support (‘‘Stiitz- 
ebenen,’’ Minkowski) sometimes called tack planes. A tack plane is a plane 
which does not separate any two points of the given point set and contains at 
least one point of this set. A plane is said to separate two points if, when the 
coordinates of the points are written in the equation of the plane two values 
with opposite signs result. These definitions enable us to find those points P, 
which lie on the boundary of (B) and to determine this boundary. (E_.g. for 
r = 3 we have to find such triples of 2, k, 1, that the determinant which represents 
the equation of the plane through these three points has the same sign for all 
possible other points P,. If the S, are the first three moments with respect to 
the origin, these determinants become Vandermond determinants and we find 
easily that the boundary planes are each passing through two neighboring 
points P,, P,4;, and one of the endpoints P; or P,. If p = 2, and the first 
two moments are given, the boundary of (B) consists of the polygon P,P; --- 
P,P, ). Then we find without difficulty the conditions to be satisfied by S in 
the form of linear inequalities between the given S;, S2, --- S,. 

We get the continuous case as a limit of the discontinuous case as t, — ty 41 
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and the points P(t) take up the whole curve (C), e.g. between t = 0 and o. 
Then the relations between the given S, become non-linear inequalities, well 
known for the problem of moments. 
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AN INEQUALITY FOR MILL’S RATIO 
By Z. W. BirNBAUM 
University of Washington 
Mr. R. D. Gordon’ recently proved the inequalities 


x 1 —}r2 1 [ —4ht2 1 —}z2 
a 6 —s e* adt< —e¢e  , forz>O0. 
z+ 1 /2n 2n vz V 2 


In the present note we show that the lower inequality can be replaced by the 
better estimate 








33 | = 








V4 + a%— 0 1 tet i in | et” dt. 
2 V/ 24 V/ 24 Jz 
Proor: According to a well-known theorem of Jensen’, for f(t) convex and 
g(t) = O in the interval (a, b), the following inequality holds 


lf tg(t) a/ | g(t) a | < [ 10 g(t) a/f g(t) dt. 


Fora = 2,b = ~, f(t) = 1/t, g(t) = te"? this inequality gives 


| war] | rewas | o a/ f te?” dt. 


Since 


z 


/ te3” dt = 6 and / fe?” dt = xe" + | e*” dt, 
z z 





1R. D. Gordon, ‘‘Values of Mill’s ratio of area to bounding ordinate of the normal 
probability integral for large values of the argument,’’ Annals of Math. Stat., Vol. 12 (1941), 
pp. 364-360. 


2 See for example: G. H. Hardy, J. E. Littlewood and G. Pélya, Inequalities, Cambridge, 
1934, p. 150-151. 
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